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Preface 


In  year  2002,  at  the  World  Conference  on  Computational  Mechanics  (WCCM)  in  Vienna,  there 
will  be  a  key  topic  for  the  first  time:  ’’Verification  and  Validation  in  Computational  Mechan¬ 
ics”  (V&V).  Besides  one  main  session  about  ’’Reliability  and  Safety”,  there  are  two  minisym¬ 
posia  about  V&V  organized.  An  accademic  V&V-minisymposium  is  organized  by  Professor 
N.-E.  Wiberg  (Europe)  and  an  industrial  one  is  organized  by  Dr.  L.  Schwer  (USA).  We  are 
participating  in  the  latter  one. 

This  final  report  -  Demonstration  of  a  Procedure  to  Establish  Guidelines  for  Reliable 
and  Safe  Numerical  Simulations  -  presented  here,  is  our  contribution  of  a  cowork  on  V&V. 
It  has  been  conducted  by  the  European  Office  of  Aerospace  Research  and  Development,  US  Air 
Force  Research  Laboratory,  US  Air  Force  Office  of  Scientific  Research. 

The  contractual  partners  are  the  University  of  the  Bundeswehr  Munich  in  Germany  -  Dept  of 
Civil  Engineering  -  Professor  N.  Gebbeken  and  the  U.S.A.  -  US  Air  Force  Research  Office.  The 
contractual  processing  time  has  been  from  June  1®*  2000  to  May  31^*  2001,  and  the  deadline  of 
this  final  report  is  July  14*^  2001. 

Our  work  has  to  be  understood  as  a  discussion  proposal.  Our  focus  of  interest  in  V&V  is 
certainly  not  terminated  by  this  work.  Review  and  comments  are  most  welcome,  and  this  report 
is  free  for  discussion  now.  Therefore,  we  are  also  looking  forward  to  an  interesting  and  lively 
discussion  at  the  WCCM  next  year. 


University  of  the  Bundeswehr  Munich 
Dept  of  Civil  Engineering  -  BauV2 
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D-85577  Neubiberg  —  Germany 
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(Professor/Dr.  Norbert  Gebbeken) 


Chapter  1 


Motivation 

1.1  Introduction 

Numerical  analysis  is  not  only  the  domain  of  experts.  Low  cost  computing  power  and  a  wide 
range  of  commercial  analysis  tools  are  readily  available.  Tools  such  as  CAD  geometry  transfer, 
automatic  meshing,  adaptive  refinement  and  optimization  devices  means  that  products  have 
become  much  easier  to  use. 

This  has  led  to  many  analyses  being  carried  out  throughout  the  entire  product  design  cycle  and 
to  the  rapid  use  of  analysis  technologies  by  smaller  companies  (NAFEMS  [27]). 

Not  only  companies,  but  also  scientific  institutions  are  tending  to  reduce  costly  and  time  con¬ 
suming  experiments  and  to  rely  more  and  more  on  the  results  of  numerical  simulations  in  order 
to  reduce  the  time  and  the  amount  of  money  needed  for  research  and  development.  Unfor¬ 
tunately,  institutions  often  have  not  the  education  needed  for  proper  use  of  these  techniques. 
Numerical  methods  are  -  methodical  inherent  -  approximate  methods,  and  therefore  they  are 
open  for  misuse  and  the  production  of  any  possible  error.  Subsequently,  model  validation  is 
necessarily  required.  Indeed,  the  cost  benefit  analysis  has  become  predominant  in  the  field  of 
research  &  development  in  engineering,  in  natural  science,  etc.  . 

There  is  a  considerable  concern  that  the  accuracy  of  numerical  simulations  using  finite  methods 
is  required  to  be  verified  in  order  to  allow  the  results  to  be  effectively  evaluated.  ’’Accuracy”, 
’’verification”,  ’’evaluation”,  (and  other  terms)  are  terms  to  be  well  defined  (section  1.3). 


1.2  Objective 

The  objective  of  this  study  is  the  demonstration  of  a  procedure  to  establish  guidelines  for 
reliable  and  safe  numerical  simulations.  Therefore,  reliable  numerical  studies  of  engineering 
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tasks  have  to  be  established  and  the  credible  use  of  finite  methods  as  well  as  the  ensurance  of 
their  accuracy  has  to  be  promoted. 

Main  fields  of  investigation 

Doing  this,  three  main  fields  of  investigation  have  to  be  addressed: 

1.  Theories  of  finite  methods, 

2.  numerical  tools  and 

3.  the  user. 

Key  questions 

Therefore,  three  key  questions  have  to  be  approached: 

1.  What  are  the  limitations  of  the  theories  of  finite  methods  used? 

2.  What  are  the  numerical  tools  able  to  solve? 

3.  What  abilities  are  required  from  the  user? 

The  main  fields  have  to  be  discussed  with  respect  to  four  levels: 

•  Evaluation, 

•  Validation, 

•  Accreditation,  and 

•  Quality  assurance, 

in  order  to  guarantee  the  quality  of  finite  methods,  software  and  user,  and  their  proper  interac¬ 
tion.  These  terms  will  be  defined  in  the  following  in  section  1.3. 

Focus 

Because  the  entire  field  of  finite  methods  and  numerical  simulations  is  so  broad  we  focus  in  this 
research  program  on  penetration  and  contact  explosion  problems. 

Nevertheless,  there  are 
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•  general  statements,  valid  in  any  case,  and 

•  specific  statements,  only  valid  for  the  problem  studied  and  the  used  software  package. 


1.3  Provision  of  a  standard 

One  intention  of  our  investigations  presented  here  is  to  provide  the  basics  for  the  reliability 
of  numerical  simulations,  in  order  to  develop  a  kind  of  standard. 

Herein,  specific  key  terms  will  be  well  defined.  Such  are,  for  example: 

•  Errors  and  user  mistakes  (section  1.3. 2. 3) 

•  Model  simplifications  (section  2.2) 

•  Code  validation  and  verification  (section  3.3) 

By  this  standard  at  hand,  the  user  will  be  able  to  evaluate  his  results,  to  detect  and  estimate 
errors  on  each  level  of  modeling  and  ultimately  to  improve  his  numerical  calculations  concerning 
their  reliability.  The  user  is  accompanied  by  a  guideline  during  all  necessary  steps  of  his 
analysis  that  is  transferable  for  his  specific  tasks. 

The  following  definitions  are  primarily  taken  from  DIN  EN  ISO  9000:2000-12  (’’quality  man¬ 
agement  systems:  fundamentals  and  vocabulary”)  [14].  Some  definitions  are  modified  here  with 
respect  to  our  specific  subject  of  interest. 

1.3.1  Scope  of  application 

Quality  management  (QM)  and  quality  assurance  are  also  relevant  here.  Specific  prin¬ 
ciples  and  key  terms  have  been  well  defined  in  the  standard  DIN  EN  ISO  9000  that  has  been 
accepted  by  the  European  committee  for  standardization  (CEN)  as  an  European  standard  in 
December  2000. 

The  basic  key  terms  will  be  defined  in  the  section  below. 

The  origins  of  the  ISO  9000  family  are  -  among  others  -  related  to  military  standards.  The 
first  ISO  (International  Standard  Organisation)  -  standards  for  quality  assurance  are  based  on 
various  national  quality  standards  of  the  80’ies  and  90’ies;  especially:  The  Canadian  standards 
”Z  299”,  the  British  standard  ”BS  5750-quality  assurance”,  the  US  military  NATO  standard 
”AQAP”  and  the  German  ”DIN”  standards. 
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DIN  EN  ISO  9000  includes  the  basics  of  QM  and  is  the  first  part  of  the  ISO  9000  family  of 
standards  that  has  been  developed  to  assist  organizations,  of  all  types  and  sizes,  to  implement 
and  operate  effective  quality  management  systems. 

This  standard  is  applicable  to  the  following: 

1.  organizations  seeking  advantage  through  the  implementation  of  a  QM  system, 

2.  organizations  seeking  confidence  from  their  suppliers  that  their  product  requirements 
will  be  satisfied, 

3.  users  of  the  products, 

4.  those  concerned  with  a  mutual  understanding  of  the  terminology  used  in  quality  manage¬ 
ment, 

5.  those  internal  or  external  to  the  organization  who  assess  the  QM  system  or  audit  it 
for  conformity  with  the  requirements  if  ISO  9001  (e.g.  auditors,  regulators,  certifica¬ 
tion/registration  bodies), 

6.  those  internal  or  external  to  the  organization  who  give  advice  or  training  on  the  QM  system 
appropriate  to  that  organization, 

7.  developers  of  related  standards. 

No.  2.,  3.  and  7.  are  especially  of  interest  here,  and  thus,  we  are  motivated  to  go  further  into 
detail. 


1.3.2  Basic  key  terms 

1.3. 2.1  General  terms  related  to  quality 

requirement 

Need  or  expectation  that  is  stated,  generally  implied  or  obligatory. 


quality 

Degree  to  which  a  set  of  inherent  characteristics  fulfills  requirements. 
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capability 

Ability  of  an  organization  or  process  to  realize  a  product  that  will  fulfill  the  requirements  for  that 
product. 


quality  assurance 

Part  of  quality  management  ( QM)  focused  on  providing  confidence  that  quality  requirements  will 
be  fulfilled. 

quality  characteristic 

Inherent  ’distinguishing  feature ’  ( —characteristic )  of  a  product  or  process  related  to  a  require¬ 
ment. 


continual  (quality)  improvement 

Increasing  the  ability  to  fulfill  quality  requirements  as  a  recurring  activity. 

effect  ivenes 

Extent  to  which  planned  activities  are  realized  and  planned  results  achieved. 

efficiency 

Relationship  between  the  result  achieved  and  the  resources  used. 


interested  party 

Person  or  group  having  an  interest  in  the  performance  or  success  of  an  organization. 

process 

Set  of  in  interrelated  or  interacting  activities  which  transforms  inputs  into  outputs. 

product 

Result  of  a  process.  There  are  four  generic  product  categories:  services,  software,  hardware, 
processed  materials. 

dependability 

Collective  term  used  to  describe  the  availability  performance  and  its  influencing  factors:  reliabil- 
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ity  performance,  maintainability  performance,  and  maintenance  support  performance.  Depend¬ 
ability  is  used  for  general  desciptions  only  and  it  is  not  quantifiable. 


traceability  and  transparency 

Ability  to  trace  the  history,  application  or  location  of  that  which  is  under  consideration. 


robustness 

The  question  of  robustness  about  a  product  is  a  criterion  for  the  sensitivity  of  the  process:  How 
do  specific  (but  minor)  modifications  within  the  process  influence  the  resulting  product?  (own 
definition). 

conformity  and  accuracy 

Conformity  is  the  fulfillment  of  a  requirement.  In  our  subject,  the  accuracy  of  numerical  results 
is  a  key  requirement  that  has  to  be  fulfilled.  (’’Non-conformity”  and  related  terms  will  be  defined 
below.) 

1.3. 2. 2  Terms  related  to  examinations 


inspection  and  testing 

Inspection  is  a  conformity  evaluation  by  observation  and  judgement  accompanied  as  appropriate 
by  measurement,  testing  or  gauging. 

verification 

Confirmation,  through  the  provision  of  objective  evidence  (data  supporting  the  existence  or  verity 
of  something),  that  specified  requirements  have  been  fulfilled. 

validation 

Confirmation,  through  the  provision  of  objective  evidence  (s.  verification),  that  the  requirements 
of  a  specified  intended  use  or  application  have  been  fulfilled. 


qualification  process 

Process  to  demonstrate  the  ability  to  fulfill  specified  requirements.  Qualification  can  concern 
persons,  products  or  systems. 
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review 

Activity  undertaken  to  determine  the  suitability,  adequacy  and  effectiveness  of  the  subject  matter 
to  achieve  established  objectives.  Example:  Nonconformity  review. 

audit 

Systematic,  independent  and  documented  process  for  obtaining  audit  evidence  and  evaluating  it 
objectively  to  determine  the  extent  to  which  (specific)  audit  criteria  are  fulfilled.  Internal  audits, 
sometimes  called  ’’first-party  audits”,  are  conducted  by,  or  on  behalf  of,  the  organization  itself 
for  internal  purposes  and  can  form  the  basis  for  an  organization  s  self-declaration  of  conformity. 
External  audits  include  what  are  generally  termed  ’’second-”  or  ’’third-party  audits  .  Second- 
party  audits  are  conducted  by  parties  having  an  interest  in  the  organization,  such  as  customers, 
or  by  other  persons  on  their  behalf.  Third-party  audits  are  conducted  by  external  independent 
organizations.  Such  organizations  provide  certification  or  registration  of  conformity  with  specific 
requirements. 

1.3.2. 3  Specific  terms  related  on  non-conformity  and  others 

At  this  point  we  would  like  to  point  out  that  we  have  to  distinguish  between  error,  mistake,  and 
other  terms  related  to  this  topic. 

Definitions: 

non-conformity 

Non-fulfillment  of  a  requirement. 

defect 

Non-fulfillment  of  a  requirement  related  to  an  intended  or  specified  use. 


preventive  action 

Action  to  eliminate  the  cause  of  a  potential  nonconformity  or  other  undesirable  potential  situa¬ 
tion. 


corrective  action 

Action  to  eliminate  the  cause  of  a  detected  nonconformity  or  other  undesirable  situation. 
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correction 

Action  to  eliminate  a  detected  non-conformity. 


error 

Error  is  a  deviation  from  the  ideal.  And  an  error  can  be  methodical  inherent.  All  finite  methods 
produce  errors.  Therefore,  a  research  field  has  been  established  that  deals  with  ’error  detection’, 
’error  estimation’  and  ’error  minimization’. 


mistake 

Mistake  is  ’’doing  wrong”.  A  mistake  is  something  that  must  be  avoided. 
Such  mistakes  are,  for  example: 


•  Wrong  use  or  a  misestimation  of  specific  methods, 


•  misapplication  of  software. 


•  misinterpretation  of  numerical  results. 


•  misinterpretation  of  errors(!). 


1.3.3  About  ’’Accreditation,  Certification  and  Quality  Assurance” 

The  following  is  an  example  of  an  general  idea.  It  is  cited  from  the  Accreditation,  Certification 
and  Quality  Assurance  Institute  [ACQuIn,  WOLPF,  Universitat  Bayreuth,  95440  Bayreuth]: 

Accreditation  is  an  evaluation  based  on  agreed  standards,  resulting  in  a  formal,  public  recognition 
of  a  programme  (or  an  institution).  It  is  a  democratic,  transparent  process  resting  upon  self-  and 
peer- assessment  for  improvement  of  academic  quality  and  public  accountability.  There  are  two 
complementary  procedures:  first  an  evaluation  procedure  consisting  of  the  self  evaluation  of  the 
applicant  plus  the  on-site  inspection  and  the  report  of  the  evaluator  group,  followed  by  the  ac¬ 
creditation  procedure  in  which  the  relevant  expert  team  analyses  and  discusses  the  self  report  and 
the  evaluators’  recommendation,  before  the  accreditation  commission  decides  on  accreditation, 
conditional  accreditation  subject  to  provisos  or  conditions,  or  denial  of  accreditation. 


CHAPTER  1.  MOTIVATION 


10 


1.3.4  Principles  of  quality  management  concerning  our  objective 

Success  can  result  from  implementing  and  maintaining  a  management  system  that  is  designed 
to  continually  improve  performance  while  adressing  the  needs  of  all  interested  parties  [14]. 
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Figure  1.1:  Interested  parties  and  quality  assurance 


As  fig.  1.1  shows,  the  interested  parties  (here)  are: 


•  method  and  code  developers 

•  tools’  vendors 

•  users  who  dispose  numerical  results 

•  internal  and  external  partners  (e.g.  testing  institutes) 

•  persons/customers  who  apply  the  results 

Of  course,  all  parties  have  specific  interests  and  requirements  that  might  differ  from  each  other. 
But  their  proper  interaction  is  crucial  for  an  efficient  quality  assurance. 

The  user,  for  example,  has  got  a  legitimate  requirement  that  the  promisses  of  the  code  vendors, 
who  primarily  want  to  sell  their  products,  are  fulfilled.  Those  who  apply  numerical  results  have 
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certain  interest  that  the  results  -  e.g.  presented  and  disposed  by  the  user  -  are  reliable  and  safe 
or  validated  and  verified.  But,  up  to  now,  the  only  one  responsible  for  the  numerical  result  is 
the  user! 

Some  of  the  quality  management  principles  identified  in  the  standard  are  basically  also  relevant 
here: 

1.  customer  focus 

2.  involvement  of  people 

3.  system  approach  to  management 

4.  continual  improvement 

5.  factual  approach  to  decision  making 

1. )  Customer  focus;  Organizations  depend  on  their  customers  and  therefore  should  under¬ 
stand  current  and  future  customer  needs,  should  meet  customer  requirements  and  strive  to  exceed 
customer  expectations. 

It  is  a  principal  part  of  QM  to  focus  on  the  customer.  Here  we  have  to  focus  on  the  per¬ 
son/organization  who  wants  to  apply  results  of  numerical  simulations  to  his/its  specific  technical 
problem.  Therefore,  the  reliablility  of  results  is  a  main  interest.  Thus,  customer  orientation  is 
also  important  here  for  all  parties  mentioned  above,  especially  for  code  vendors. 

Considered  the  code  vendor  and  the  code  user  as  contractual  partners,  customer  requirements 
have  to  be  defined  and  satisfied  (fig.  1.1).  The  user  requires  tools  with  characteristics  that 
satisfy  his  needs  and  expectations.  This  needs  and  expectations  are  expressed  in  product  speci¬ 
fications  and  collectively  referred  to  as  customer  requirements.  ’’The  reliable  and  safe  numerical 
simulation  of  high  velocity  impact  and  perforation  in  the  range  of  100  -  1000  ^  impact  veloc¬ 
ity,  1-100  g  (metal-)  penetrator  weight  into  a  (unreinforced)  mortar  kind  of  material”  may  be 
specified  conctractually  by  the  customer  (=  user)  or  may  be  determined  by  the  organization 
(tool’s  vendor)  itself.  In  either  case,  the  customer  ultimately  determines  the  acceptability  of  the 
product. 

Annotation:  Considered  the  code  user  and  the  persons  who  apply  results  as  contractual 

partners,  this  concept  can  be  transfered  respectively. 

2. )  Involvement  of  people:  People  at  all  levels  are  the  essence  of  an  organization  and  their 
full  involvement  enables  their  abilities  to  be  used  for  the  organization’s  benefit. 
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This  principle  can  be  transfered  to  our  task  when  education,  experience  and  competency  is  taken 
into  account  of  all  people  involved.  In  this  context,  special  focus  on  user  requirements  has  to  be 
done.  But  also  the  proper  interaction  of  people  (experts,  developers,  users,  ...)  comes  within  the 
limits  of  this  topic:  Safe  and  relible  numerical  simulations  are  -  at  least  -  the  result  of  teamwork 
based  on  the  proper  involvement  of  all  parties. 

3. )  System  approach  to  management:  Identifying,  understanding  and  managing  inter¬ 
related  processes  as  a  system  contributes  to  the  organization’s  effectiveness  and  efficiency  in 
achieving  its  objectives. 

Transfered  to  our  task,  the  intention  of  this  principle  is  to  understand  the  numerical  simulation 
as  a  process.  The  desired  result  (reliable  numerical  simulations)  is  achieved  more  efficiently 
when  activities  (method  development,  programming,  implementation  of  numerical  features  and 
models  into  the  code,  numerical  setup  and  calculation,  ...,  interpretation  of  results,  technical  ap¬ 
plication,  research  and  further  development,  ...)  and  related  recources  (monetary  environment 
and  input,  development  time,  deployment  and  personnel  qualification  v-  user  requirements,  ...) 
are  managed  as  a  process  within  a  system.  Therefore  the  partners  have  not  only  to  be  treated 
as  contractual  partners,  as  described  in  no.#l  and  the  term  ’’organization”  may  be  ’ascended’  to 
a  higher  level:  It  should  include  the  whole  system  of  all  parties  and  their  (proper!?)  interactions. 

4. )  Continual  improvement:  Continual  improvement  of  the  organization’s  overall  perfor¬ 
mance  should  be  a  permanent  objective  of  the  organization. 

Beyond  the  aim  to  improove  method,  tools  and  users  a  QM  system  has  got  the  aim  of  continual 
improvement  of  itself  by  increasing  the  probability  of  enhancing  the  satisfaction  of  customers 
and  other  interested  parties.  Analysing  and  evaluating  the  existing  situation  is  a  first  action  for 
improvement  that  has  to  be  transfered  also  to  our  topic. 

A  goal  of  QM  here,  is  the  continual  improvement  of  performance  of  methods,  codes  and  tools. 
Therefore,  the  responsible  persons  of  this  three  key  levels  have  to  separately  review  and  as¬ 
sure  quality  continuously.  The  method  has  to  be  scrutinized,  codes  have  to  be  verified  and 
user  requirements  have  to  be  defined  and  accredited.  Audits  on  each  level  have  to  be  done 
independently. 

Feedback  from  ’customers’  (see  above)  and  other  interested  parties,  audits  and  review  of  the 
QM  system  can  also  be  used  to  identify  opportunities  for  (continual)  improvement. 
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5.)  Factual  approach  to  decision  making:  Effective  decisions  are  based  on  the  analysis  of 
data  and  information. 

To  heed  this  principle  is  important  to  act  on  the  maxim  of  quality  assurance:  The  factual 
approach  must  be  fundament  of  the  evaluation  of  methods,  codes  and  users. 

In  this  context,  the  importance  of  a  documentation  of  the  entire  process  has  to  be  mentioned: 
Documentation  enables  communication  of  intent  and  consistence  of  action: 

•  achievement  of  conformity  to  stated  requirements, 

•  provision  of  appropriate  user  training, 

•  repeatability,  traceability  and  objective  evidence  of  numerical  results, 

•  evaluation  of  the  effectiveness  and  continuing  suitability  of  the  QM  system. 

For  the  purpose  of  this  definition,  documentation  is  a  value- adding  activity  in  terms  of  a  continual 
quality  improvement. 

Again,  all  levels  have  to  be  embraced:  A  documentation  of  the  method  yields  insight  in  the 
theoretical  background.  The  documentation  of  codes  prevents  ’black  boxes’.  The  documentation 
of  numerical  studies  detailed  by  the  user  is  an  instrument  of  self  control  and  contributes  to 
repeatability  of  results. 

And,  besides  the  regard  to  each  level  for  its  own,  the  association  of  the  levels  is  important,  too: 
Therefore,  essentially  the  following  types  of  documents  are  used  in  QM  systems: 

•  quality  manuals  that  provide  consistent  information  about  the  organization’s  QM  system; 

•  quality  plans  that  describe  how  the  QM  system  is  applied  to  a  specific  product  (code); 

•  specifications  that  state  requirements; 

•  guidelines  or  work  instructions  that  state  recommendations  or  suggestions  or  that  provide 
information  about  how  to  perform  activities  and  processes  consistently; 

•  (and  others) 

In  general:  The  extent  of  documentation  depends  on  factors  such  as  the  type  and  size  of 

organization,  the  complexity  and  interaction  of  processes,  the  complexity  of  products,  customer 
requirements,  the  applicable  regulatory  requirements,  the  demonstrated  ability  of  personnel. 


CHAPTER  1.  MOTIVATION 


14 


and  the  extent  to  which  it  is  necessary  to  demonstrate  fulfillment  of  quality  management  system 
requirements. 

Transfered  to  our  specific  topic,  main  focus  should  be  done  on  the  regard  to  each  level  for  its 
own.  We  are  going  further  into  detail;  see  below. 

1.3.5  Conclusion 
Generally: 

Section  1.3  has  outlined  quality  aspects  and  quality  key  terms  (section  1.3.2),  generally  and  with 
special  focus  on  our  subject. 

Without  going  deeper  into  the  details  of  standards  like  the  mentioned  DIN  EN  ISO  9000  family, 
we  asserted  that  a  proper  quality  management  is  also  necessary  for  reliability  and  safety  of 
numerical  simulations.  Again,  we  would  like  to  point  out  the  importance  of  a  proper  interaction 
of  all  interested  parties  (section  1.3.4,  fig.  1.1). 

Our  special  focus  here: 

This  report  establishes  a  guideline  for  the  user  to  detect  and  evaluate  errors  on  each  level  of  the 
numerical  analysis. 

— )•  Errors  must  be  detected  and  minimized. 

The  proposal,  presented  here,  might  also  to  be  understood  as  a  trainee  program  or  a  qualifica¬ 
tion  process  related  to  quality  assurance  for  numerical  calculations  with  special  respect  to 
hyper  dynamics. 

— >  Mistakes  will  be  reduced  by  following  these  proposals. 

The  user  is  educated  by  that  way  and  might  be  accredited  as  an  expert  of  specific  numerical 
methods.  Therefore  a  test  -  possibly  as  a  part  of  internal  or/and  external  audits  -  might  be 
included  in  the  future. 

Risks  are  unfolded  and  a  deep  insight  into  the  method  or  the  code  used  is  established,  and  this 
allows  to  see  behind  the  curtain.  The  intention  is  a  continual  improvement  of  methods,  codes 
and  users  regarding  their  quality. 
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1 . 4  Classifications 

The  following  is  mainly  related  to  the  specific  problem  of  penetration  and  contact  detonation. 
In  order  to  select  the  proper  methodology  we  need  to  define  and  classify  the  problem  precisely 
with  respect  to  load,  structure,  material  and  the  goal  of  investigation. 

1.4.1  Static  case,  structural  dynamics  and  hyper  dynamics 

First  of  all,  we  need  to  distinguish  between  time  effects  where  mass  forces  are  generated  (dynamic 
case)  and  where  mass  forces  are  not  generated  (static  case).  Furthermore,  the  dynamic  case  has 
to  be  subdivided  into  hyper  dynamics  and  ordinary  dynamics  -  i.e.  structural  dynamics.  This 
has  to  be  taken  into  account  when  using  numerical  methods: 

•  Static  case:  Analytical  or  ’hand’-  calculations  are  in  some  cases  possible  and  often  suffi¬ 
cient  to  solve  simple  static  problems  or  they  can  (roughly)  verify  numerical  calculations; 
linear /nonlinear  -  ordinary  -  FEM  calculations  are  still  necessary  and  challenging  for  large 
and/or  irregular  geometric  problems;  material  models  are  known  for  most  instances;  ref¬ 
erence  tests  and  benchmarks  are  broadly  available.  The  numerical  methods  and  the  codes 
used  are  overall  known. 

•  Structural  dynamics:  A  reference  to  the  static  case  is  oftentimes  possible  (e.g.  by 
enhancement  factors);  which  finite  method  is  best,  depends  on  the  task  to  be  solved  (see 
section  2.3.5).  The  numerical  methods  and  the  codes  used  are  state  of  the  art. 

•  Hyper  dynamics:  Analytical  calculations  are  not  available  and  ’hand’  formulas  -  if 
available  -  are  normally  confined  to  specific  empirical  studies;  to  obtain  material  data  is 
often  a  difficult  task.  The  numerical  methods  are  specialized  (e.g.  wave  propagation  codes 
/  hydrocodes)  and  some  further  error  possibilities  are  -  methodical  inherent  -  introduced 
(as  shown  below  in  chapter  3).  Therefore,  primarily  in  hyper  dynamics,  code  validation  k 
verification,  reference  tests  and  benchmarks  are  principal  duties  in  order  to  obtain  reliable 
numerical  results. 

For  the  sake  of  a  classification  of  dynamic  and  hyper  dynamic  effects  on  structures,  different  load 
types  and  important  values  as  impact  duration,  impact  velocity,  strain  rate,  ratio  of  dynamic  to 
static  strength,  material  behavior  and  related  testing  device  are  compared  in  table  1.1. 

Table  1.1  is  based  on  data  from  different  publications  that  have  been  collected  and  complemented 
here  (e.g.  Bischoff  k  Perry  (1995)  [6],  Zukas  (1982)  [36]). 
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Table  1.1:  Comparison  of  important  key  values,  taken  (and  modified)  from  Gebbeken  &:  Ruppert  [21] 
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Loadings  are  mostly  not  of  static  nature.  In  fig.  1.2,  different  load  types  and  time  dependent 
phenomena  are  related  to  measured  strain  rates  in  concrete.  In  addition,  the  increase  of  the 
compressive  strength  with  increasing  strain  rates  is  shown  up  to  the  limit  of  10^^.  There  is  still 
a  lack  of  data  beyond  that  limit. 

Not  mentioning  when  strain  rates  are  greater  than  10^^  is  not  reliable;  it  is  not  an  error,  it  is  a 
mistake,  or  at  least  ’wantonly  negligent’. 

In  fig.  1.3,  tension  stress  vs.  strain  is  shown  for  two  different  strain  rates.  These  curves  are  valid 
for  reinforcing  steel  BSt420/500,  DIN1045  (German  standard  for  reinforced  concrete  structures), 
with  a  deterministic  yield  limit  of  ayieia  =  420MPa.  The  higher  the  rate,  the  more  stress  can  be 
taken. 

Figures  1.2  and  1.3  are  two  examples  that  show  typical  material  behavior  due  to  different  loading 
situations  and  strain  rates.  What  is  shown  here  for  concrete  and  a  specific  reinforcing  steel  from 
DIN1045,  is  generally  valid  for  other  materials,  too. 

For  highly  dynamic  loads  we  also  refer  to  G.  R.  Johnson’s  publication  [25]  where  data  for 
metals  up  to  e  =  10^^  are  presented. 
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CEB  Recommendations  (20  MPa)  «  / 


Figure  1.2:  Compressive  strength  versus  strain  rate,  /c(e),  and  related  phenomena  (based  on 
Bischoff  &  Perry  [6]) 

The  behavior  of  materials  and  structures  subjected  to  static  or  quasistatic  loads,  as  well  as 
creep  &:  shrinkage  phenomena  have  been  investigated  in  the  past  decades.  Structural  dynamics 
(e  <10^^)  is  well  known  from  the  methodical  point  of  view,  but  is  still  challenging  when  applying 
new  materials  or  damping  devices,  respectively.  There  is  a  lot  of  literature  available  to  solve 
such  conventional  problems. 

Most  of  the  material  parameters  are  actually  rate  dependent.  There  are  two  ways  to  take  rate 
dependency  into  account: 

1.  rheological  modeling, 

2.  multiply  the  static  values  by  an  amplification  factor. 

The  latter  is  a  common  strategy  and  will  also  be  applied  here.  The  amplification  factor  fc^^/ fc 
(Table  1.1)  can  be  obtained  only  by  testing.  We  have  taken  the  values  of  the  amplification  factor 
from  the  trilinear  regression  curve  suggested  by  the  CEB-Bulletin  No. 187(1987),  [8]  -  [11]  as 
shown  for  concrete  in  fig.  1.2. 
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Conclusion 

Different  loading  situations  (tab.  1.1)  cause  different  material  and  structural  behavior  that  re¬ 
quire  different  numerical  methods  to  be  used. 

An  extrapolation  from  structural  dynamics  to  hyper  dynamics  is  often  not  the  appropriate  way 
(fig.  1.2),  because  of  both  lack  of  data  and  different  material  and  structural  behavior. 

In  the  following  section,  we  focus  on  hyper  dynamics. 


1.4.2  Levels  of  material  modeling  in  continuum  mechanics 

Before  going  into  details  of  material  modeling  in  hyper  dynamics,  one  should  be  aware  of  the 
different  levels  of  material  models  that  might  be  possible. 

It  depends  on  the  problem  to  be  solved,  which  level  is  the  most  appropriate  one.  In  order 
to  model  real  structures,  often  the  homogenization  hypothesis  is  adopted  which  enables  the 
formulation  of  the  constitutive  equations  on  the  level  of  macromechanics  (fig.  1.4). 
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real  microscopic  material 


simplified  microstructure 

i  i  I  i 


material  formulation 
on  the  level  of 
micromechanics 


homogeneous  microstructure 

-‘V 

classical  continuum  mechanics 


microplane 

formulations 


macroscopic  models 

-  theory  of  elasticity 

-  theory  of  plasticity 

-  viscoelasticity  -plasticity 

-  damage  and  fracture 


application  to  structural  analysis 


Figure  1.4:  Homogenization  of  a  microscopic  material  to  a  macroscopic  scale  (Ruppert  Sz  Gebbeken 
[31],  [22]) 

Annotation 

Adopting  the  homogenization  assumption  we  must  realize,  that  the  information  of  the  indivi¬ 
dual  constituents  is  smeared  out  (homogenised) .  Although  concrete,  for  example,  is  a  highly 
inhomogenious  material,  this  is  commonly  done  when  modeling  on  the  macroscopic  scale. 

Micromechanic  approximations  are  very  time  consuming,  what  sometimes  is  unacceptable;  espe¬ 
cially  for  huge  problem  setups. 

1.4.3  Structural  behavior 

Structural  behavior  is  influenced  by  many  factors  including  the  distribution  of  loads,  the  duration 
at  which  loads  are  applied,  interaction  with  adjacent  structural  components,  the  material  prop¬ 
erties,  the  stiffness  of  individual  members  and  devices  and  stiffness  of  overall  structure,  mass 
distribution,  temperature  influence,  etc.  .  All  of  these  aspects  comprise  specific  idealizations 
(simplifications)  to  some  ideal  form  representing  the  reality. 

Whether  the  impact  causes  mostly  local  (fig.  1.5(a))  and/or  global  eff'ects  (fig.  1.5(b,c)),  is 


CHAPTER  1.  MOTIVATION 


20 


essentially  based  on  peak  pressure,  impact  duration  and  area,  and  the  proportion  between  the 
masses  of  impactor  and  target. 


Figure  1.5:  Frame,  structural  response:  a)  hyper  dynamics,  b)  dynamics,  c)  quasi  statics 

An  impact  duration  within  milli-  or  microseconds  is  related  to  hyper  dynamics.  The  global 
structural  response  is  within  some  Hertz;  this  means  some  oscillations  per  second.  Therefore,  it 
becomes  evident  that  the  structure  is  often  too  inert  to  react  in  total  as  a  whole  structure  under 
hyper  dynamic  actions  (fig.  1.5(a)).  This  is  demonstrated  by  an  impactor  with  a  relative  small 
mass  penetrating  the  structure  by  high  velocity  impact.  The  heavy  structure  is  much  too  inert 
to  react  in  total.  No  observable  displacement  occurs  (w=0)  and  almost  no  response  forces  are 
generated  (N<<F/2). 

In  case  b)  and  c),  the  impact  velocity  is  less  than  in  case  a)  and/or  the  impactor  is  heavier. 
Thus,  in  addition  to  local  damage  also  global  structural  response  occurs  and  reaction  forces  in  the 
system  are  generated.  In  case  a)  and  b)  special  finite  methods  have  to  be  used,  e.g.  hydrocodes 
(chapter  3).  In  case  c)  a  classical- method-based  calculation  in  the  range  of  structural  dynamics 
(e.g.  Biggs  [5],  Eibl  [16])  may  be  sufficient. 

1.4.4  Consequences 

Section  1.4  has  shown  significant  variations  in  material  as  well  as  in  structural  responses  caused 
by  different  load  situations.  Different  engineering  tasks  ask  for  different  numerical  methods  and 
different  software  tools.  In  the  next  section,  specific  methods  are  generally  and  briefly  presented. 

The  more  challenging  the  task  investigated,  the  more  fundamental  knowledge  of  physics  and 
methods  used  is  required.  Therefore  an  introduction  into  finite  methods  will  be  presented. 
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Finite  methods  -  Relevant  basics 

2.1  Introduction 

Nowadays,  one  cannot  imagine  any  numerical  calculation  in  the  fields  of  applied  mechanics 
without  using  finite  methods.  Finite  methods  (FM)  means  finite  element  methods  (FEM)  and 
finite  difference  methods  (FDM).  The  boundary  element  method  (BEM)  and  the  family  of 
meshless  methods  (Smoothed  Particle  Hydrodynamics,  Element-Free  Galerkin,  ...)  will  not  be 
discussed  here. 

Numerical  methods  allow  us  to  replace  a  set  of  differential  equations  by  a  set  of  linearized 
equations.  The  finite  methods  are  generally  procedures  for  approximating  the  set  of  differential 
equations.  Numerical  methods  are  used  because  they  allow  approximate  solutions  of  a  broad 
range  of  problems  which  otherwise  could  not  be  solved.  Problems  with  irregular  boundaries, 
sophisticated  material  models,  geometrical  nonlinearities  and  difficult  boundary  conditions,  for 
example,  are  unsolvable  analytically.  Furthermore,  optimization  strategies,  extensive  parameter 
studies,  etc.  are  impossible  without  numerical  studies. 

The  growing  use  of  powerful  hardware  enables  the  implementation  of  more  complicated  and 
larger  systems  into  the  numerical  code  (software).  As  computing  power  became  more  widely 
available,  industry  increasingly  started  to  solve  practical  engineering  problems  using  finite  meth¬ 
ods  in  order  to  ensure  quality  and  to  reduce  time  needed  for  development. 


2.2  About  Modeling 

2.2.1  Basic  remarks 

Before  one  executes  a  numerical  simulation  by  using  finite  method  codes,  a  mechanical  model 
which  represents  the  reality  has  to  be  implemented  into  the  code  by  numerical  modeling. 
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In  general,  the  different  steps  of  modeling  are  roughly  as  follows; 

Reality  — )■  simplifications  — >  Mechanical  model  -A-  simplifications  -4  Numerical  model. 

Although  computer  resources  now  permit  models  of  high  complexity  to  be  processed,  in  most 
analyses  there  remains  a  requirement  to  simplify  the  problem  to  be  solved.  Simplifications  are 
necessary  ’’engineering  approaches”. 

Every  model  is  an  abstraction  and  simplification  of  the  reality,  and,  therefore,  incorporates 
model-inherent  errors;  deviations  from  the  reality. 

The  degree  of  simplification  depends  on  both,  the  accuracy  required  and  the  resources  available. 
Thus,  when  the  idealized  structure  is  analysed,  by  finite  element  methods  for  example,  deviations 
from  the  real  structure  will  certainly  arise  depending  on  the  assumptions  made. 

The  following  basic  remarks  are  taken  (and  modified)  from  Reddy  et  al.  [29]. 

Modeling  is  an  art,  based  on  the  ability  to  ’’visualize”  the  physical  reality.  All  basic  and  applied 
knowledge  of  physical  problems,  finite  elements,  and  solution  algorithms  contributes  to  modeling 
expertise. 

In  (numerical)  modeling,  the  principal  difficulty  faced  by  a  typical  user  of  a  code  is  to  understand 
the  physical  action  and  boundary  conditions  of  the  actual  structure,  and  the  limitations  of  ap¬ 
plicable  theory,  well  enough  to  prepare  a  satisfactory  model.  Another  difficulty  is  to  understand 
the  behaviors  of  various  elements,  and  the  tool’s  options  or  code’s  limitations,  well  enough  to 
make  an  intelligent  choice  among  them.  The  result  may  be  a  poor  specification  of  the  problem 
to  be  solved,  a  model  that  fails  to  reflect  important  features  of  the  physical  problem,  fine  detail 
irrelevant  to  the  problem,  a  solution  based  on  inappropriate  loading  or  support  conditions,  and  a 
surplus  of  computed  results  which  are  not  properly  examined  and  questioned  (I.C.  Taig:  ’’Finite 
Element  Analysis  in  Industry  -  Expertise  or  Profiency?”,  1984). 

Numerical  modeling  is  more  than  just  laying  out  a  mesh. 

2.2.2  The  way  in  general 

Fig.  2.1  -  taken  (and  modified)  from  K.-J.  Bathe  ’’Finite-Elemente-Methoden”  [2]  -  summarizes 
the  process  of  finite  method  analysis  and  it  shows  the  way  in  general. 

The  finite  method  will  solve  only  the  model  selected.  All  assumptions  made  for  this  model  will 
be  reflected  therefore.  Hence,  the  choice  of  an  appropriate  mechanical  model  -  or  some  specific 
models  -  is  a  key  step,  and  the  solution  of  the  mechanical  model  is  only  one  part  of  the  process 
of  analysis. 
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Figure  2.1:  The  way  of  finite  method  analysis  in  general;  taken  (and  modified)  from  Bathe  [2] 


2.2.3  Introductory  examples 


The  way  described  before  was  a  general  way  of  the  process  of  analysis.  Nevertheless,  one  should 
carefully  specify  the  scheme  shown  in  fig.  2.1.  The  authors,  e.g.,  have  specified  the  scheme  with 
respect  to  contact  detonation  problems  in  [32],  that  is  abstracted  in  section  2. 2. 3. 3. 
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2. 2. 3.1  Mechcinical  modeling  -  example:  Structural  engineering  -  steel  frame 


First  of  all,  some  pre-work  has  to  be  done  in  order  to  find  the  proper  mechanical  model.  Such  can 
be,  for  example,  hand  calculations,  previous  analysis  on  smaller  and  easier  models  and  physical 
tests.  Doing  so,  a  careful  appraisal  of  FE  results  gives  more  insight  into  model  behavior. 

A  physical  task  might  be,  for  example,  to  design  a  steel  frame  construction.  The  actual  structure 
is  shown  in  fig.  2.2.  Fig.  2.3  illustrates  specific  mechanical  models  of  the  given  problem  of  fig.  2.2 
as  given  on  a  shop- drawing. 

Indeed,  this  task  is  understood  as  a  fundamental  engineering  job,  but  already  up  to  this  point 
inevitable  errors  (here:  deviations  from  reality)  occur,  because  an  exact  and  unequivocal  de¬ 
scription  of  reality  does  not  exist.  That  is  why  fig.  2.3  illustrates  only  a  choice  of  some  feasible 
mechanical  idealizations  that  are  typically  used  in  the  field  of  engineering  mechanics  and  struc¬ 
tural  engineering.  A  basic  lecture  at  the  university  on  structural  engineering  might  contain  these 
illustrations. 


Figure  2.2:  Example:  Real  frame  structure 


O 


Figure  2.3:  Example:  Idealization  of  frame 
structure  -  specific  mechanical  models 


Note:  We  have  chosen,  in  a  straight  forward  manner,  beam  elements.  Why  not  shell  el¬ 

ements;  what  is  obout  the  structural  connectors,  the  fasteners,  the  weldings,  multidirectional 
(3d-)  influences,  etc.  etc.  ? 

After  a  useful  mechanical  model  has  been  selected,  the  next  step  is  to  find  the  proper  loading 
cases  and  conduct  a  design  procedure  in  accordance  with  the  specific  national  standards.  As  well 
as  the  mechanical  model  contains  idealizations  there  are  also  approximations  and  engineering 
assumptions  within  the  loads,  foundations,  solid  conditions,  structural  connections,  assembling 
process,  etc.  etc.  . 
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2. 2. 3. 2  About  engineering  design 

It  is  worth  mentioning  that  there  is  a  difference  between  engineering  design  tasks  (e.g.  the 
problem  discussed  above)  and  reliability  tasks  where  the  goal  of  investigation  is  a  research 
project  that  is  targeted  on  exactness. 

Modeling  is  also  part  of  the  design  process.  Engineering  design  is  related  to  statistical  consider¬ 
ations  concerning  the  safety  of  structures.  Design  engineers  have  to  make  sure  that  the  ultimate 
limit  of  load  capacity  (and  also  the  impairment  of  serviceability  and  durability)  of  the  structure 
is  not  achieved.  This  is  ensured  by  specific  safety  factors  in  the  (national)  standards. 

Simplifications,  e.g.  in  the  mechanical  model,  may  be  accepted  during  the  engineering  design 
process,  if  we  are  on  the  safe  side  and  the  structure  is  still  economical  (or:  Safety  is  payed!). 

But  also  here,  the  question  is:  How  do  model-decisions,  restrictions  and  simplifications  influence 
the  level  of  safety  according  to  the  postulated  safety  of  the  standards?  The  competition  of 
engineering  companies  often  drive  us  to  approach  the  postulated  safety  of  standards  as  best  as 
possible.  Large  safety  reserves  are  too  expensive. 

And  therefore,  it  is  important  to  ask  how  to  get  safe  and  reliable  numerical  simulations  based 
on  mechanical  models  that  are  close  to  reality  at  as  best  as  possible. 

Some  examples  of  wrong  design 

Errors  in  engineering  design  can  occur.  Design  by  using  finite  element  analysis  is  also  prone 
to  error.  In  [’’Application  of  Computer  Technology”  by  P.  Gardener,  Forensic  Engineering 
Conference,  1999.]  a  survey  carried  out  by  the  American  Society  of  Civil  Engineers  (ASCE) 
is  mentioned.  The  survey  covered  52  cases  of  error  where  the  design  process  involved  use  of 
computers.  It  was  found  that  17%  of  the  cases  were  affected  by  hardware  error,  25%  by  software 
error  and  the  remaining  58%  by  some  form  of  user  mistake: 

1.  hardware  errors  (17%), 

2.  software  errors  (25%), 

3.  user  mistakes  (58%). 

’’Some  form  of  user  mistakes”  can  be  itemized  as  follows: 

•  wrong  mechanical  models  implemented  by  the  user. 
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•  the  ’’wrong  method  on  the  wrong  place”, 

•  ’’simple  input  errors”:  input  mistakes,  oversights,  etc., 

•  last  but  not  least:  misinterpretations  of  the  numerical  output. 

These  user  mistakes  are  on  the  one  hand  mainly  based  on  the  lack  of  fundamental  knowledge, 
and  on  the  other  hand  based  on  time  pressure  and  rush. 


2. 2. 3. 3  Numerical  modeling  -  introductory  example  for  hyperdynamic  tasks 


Special  focus  is  done  on  hyper  dynamics.  An  example,  therefore,  is  a  contact  detonation  problem 
that  will  be  described  in  detail;  see  section  3.5.1. 


Figure  2.4:  Modeling  of  a  contact  detonation  problem,  [32] 


The  procedure  of  modeling  is  shown  in  figure  2.4  (taken  and  modified  from  [32]). 
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In  chapter  4  a  more  detailed  guideline  will  be  presented  that  accompanies  the  user  for  -  amongst 
other  things  -  the  sake  of  proper  modeling. 


2. 2. 3. 4  An  input  check  list  (in  general) 

The  following  list  presents  a  general  FE  analysis  input  check  list,  taken  (and  modified)  from 
[BENCHmark  journal  of  NAFEMS,  07/2000].  It  was  not  the  goal  here,  to  go  into  detail  at  every 
stage  of  the  design  process. 


Table  2.1:  input  check  list  for  a  general  FE  ananysis 


analysis  type 

elements  missing 

units 

elements  duplicated 

extent  of  model 

consistent  normals 

material  data 

constraint  equations 

co-ordinate  system 

symmetry  constraints 

major  dimensions 

supports 

element  types 

boundary  conditions 

real  constants 

load  cases 

mesh  density 

summed  mass 

element  plots 

master  freedoms 

element  shapes 

front/band  width 

internal  edges 

output  options 

The  check  list  in  tab.  4.3  is  general  and  not  precise  enough  for  specialized  problems.  But  it  is  a 
guidance  and  it  has  to  be  modified  and  supplemented  with  respect  to  the  problem  that  has  to 
be  solved  (chapter  4). 


2.2.4  Conclusion 

Section  2.2  has  -  among  other  things  -  indicated  the  sources  of  errors  that  might  be  caused 
during  the  modeling  process.  In  the  following  special  focus  is  given  to  sources  of  errors  within 
the  method  (section  2.3)  and  to  simplifications  concerning  material  modeling  (section  2.4). 
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2.3  Sources  of  errors  within  the  method 

2.3.1  Finite  methods  -  Introductory 

Finite  methods  are  approximation  methods,  because  a  variational  principle  adopting  trial  func¬ 
tions  is  used  yielding  a  set  of  equations,  here  the  equilibrium  formulation  in  FEM: 

K(v)-v-F  =  '_!0^  (2.1) 

The  stiffness  matrix  K(^)  is  assembled  representing  both,  material  and  system  properties;  v 
represents  the  unknown  vector  of  nodal  displacements;  the  load  vector  F  has  to  be  given  by  a 
specific  load  condition. 

The  nil  ("0")  in  eq.  2.1  represents  a  mathematical  truncation.  It  is  dependent  on  the  defined 
criteria  of  how  becomes  almost  0. 


Figure  2.5:  Load-displacement-path,  internal  truncation  algorithm  of  FEM 

Fig.  2.5  shows  the  trace  of  an  incremental  iterative  algorithm  driven  by  adaptive  increments 
calculated  from  defined  ’’norms”.  This  ’’internal  truncation  algorithm”  is  based  on  the  principle 
of  exactness  (section  3.4.3). 

In  the  following,  specific  groups  are  listed  that  -  at  each  case  -  represent  error  levels: 

•  structural  component  (rods,  plates,  ...), 

•  geometry  of  the  numerical  idealization, 


•  material  formulation 
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•  contact  algorithms, 

•  description  of  frictional  effects, 

•  treatment  of  distorted  elements, 

•  specific  items,  e.g.  of  hyper  dynamics: 

—  shock  oscillations, 

—  Lagrange/Lagrange  interaction, 

—  Euler/Lagrange  coupling, 

—  treating  the  multimaterial  processor, 

—  time  step  regulation  in  the  explicit  integration  scheme. 

On  each  level,  possible  errors  have  to  be  detected,  treated  and  quantified  sepa¬ 
rately! 

2.3.2  A  general  example  for  approximation  methods 

As  a  first  general  example  for  approximation  methods,  the  calculation  of  the  circumference  of  a 
circle  will  shown  {taken  (and  modified)  from  J.N.  Reddy  et  al.,  [30]}: 


The  problem  to  estimate  the  perimeter  of  a  circle  of  radius  r  might  be  solved  by  approximating 
it  by  line  segments,  whose  lengths  have  to  be  summarized.  This  is  a  trivial  example,  but  it 
illustrates  some  crucial  ideas  and  steps  involved  in  a  finite  element  analysis. 

The  first  step  is  the  ’’discretization”:  The  domain  is  the  circumference  of  a  circle  that  is  rep¬ 
resented  as  a  collection  of  a  finite  number  n  of  subdomains.  These  n  line  segments  are  called 
’’elements”  and  the  collection  of  them  is  called  ’’finite  element  mesh”.  The  elements  are  con¬ 
nected  to  each  other  at  their  points;  so  to  call  ’’nodes”.  A  special  case  of  this  discretization  is 
the  uniform  one,  when  all  segments  are  of  the  same  length. 

The  second  step  is  to  formulate  the  ’’element  equations”.  Therefore  a  typical  element,  i.e.  line 
segment,  is  isolated  and  its  required  properties,  i.e.  length,  are  computed  by  some  appropriate 
means.  Let  he  be  the  length  of  element  in  the  mesh: 

he  =  2r  ■  sin^cpe  ,  (2'2) 

where  r  is  the  radius  of  the  circle  and  (j)e  <  tt  is  the  angle  subtended  by  the  line  segment. 
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In  the  third  step  an  assembly  of  element  equations  and  solutions  have  to  be  done.  The  approxi¬ 
mate  value  of  the  circumference  (or  perimeter)  of  the  circle  is  obtained  by  putting  together  the 
element  properties  in  a  meaningful  way;  this  process  is  called  the  ’’assembly”  of  the  element 
equations.  In  our  example,  it  is  computed  by  the  sum  of  the  lengths  of  all  elements: 

f/n-Ehe.  (2.3) 

e=l 

Then  Un  represents  an  approximation  to  the  actual  perimeter  u  by  using  n  segments. 


If  the  mesh  is  uniform,  (j)e  =  ^  and  he  is  always  the  same  for  each  of  the  segments  in  the  mesh. 
Then  Un  becomes: 


(2.4) 


For  this  simple  problem,  we  know  the  exact  solution: 


u  =  2'Kr  . 


(2.5) 


Thus,  the  ’’error”  En  in  the  approximation  can  be  estimated  by  the  following  term: 

(Stt  tt  'N 

- 2sin—  I  =  2nr  —  Un  ■ 

n  nj 


(2.6) 


Here,  En  goes  to  zero  as  n  — >■  oo  that  can  be  shown  as  follows.  Letting  a:  =  ^,  we  have 


Un 


^  TT  „  simrx 

2r  ■  n  ■  sin—  =  2r - 

n  X 


and 


,  simrx\  U  cos'!rx\ 

hm  Un  =  hm  2r - )  =  hm  2'Kr -  =  27rr  =  u  . 

x^oo  x-^0  \  X  /  x->0  \  1  / 


[30] 


Annotation: 

This  calculation  is  well  understood  as  an  ’’asymptotic  convergence  study”.  For  this  specific 
example,  the  circumference  of  a  circle  can  be  approximated  as  closely  as  we  wish  by  a  finite 
number  of  piecewise-linear  functions.  As  the  number  of  elements  (segments),  is  increased,  the 
approximation  improves  and  the  error  in  the  approximation  decreases  (fig.  2.6,  E— >0). 

Certainly  this  circumference  approximation  is  very  simple,  but  this  presentation  has  revealed 
some  specific  problems  and  steps  that  are  also  important  when  talking  about  finite  methods. 
The  additional  problems  are  described  in  detail  in  the  following  chapters  depending  on  the 
method  used. 
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Figure  2.6:  Asymptotic  convergence  of  a  simple  circumference  approximation 

1.  It  is  up  to  the  user  to  define  the  Error  E  which  can  be  accepted.  From  that  depends  the 
’’size”  of  the  problem. 

2.  If  an  element  is  too  stiff,  too  weak,  or  shows  locking,  the  solution  might  asymptotically 
converge  to  a  wrong  solution! 

What  is  the  proper  method? 

Uncounted  papers  and  textbooks  have  been  published  dealing  with  theory  and  application  of 
FE  or  FD  methods.  There  are  several  computer  codes  which  are  commercially  available  to 
run  on  wide  range  of  computer  hardware.  Meanwhile,  FE  analysis  is  playing  a  crucial  role  in 
research,  development  and  practice.  Computer  simulation  has  become  an  important  phase  of 
design  processes. 

The  key  question  is:  ’’What  is  the  proper  method?”.  That  leads  to  a  comparison  of  different 
methods  in  order  to  gain  information  about  both,  advantages  and  disadvantages. 

2.3.3  The  finite  element  method  (FEM) 

2.3.3. 1  Visible  description 

To  apply  the  finite  element  method,  a  region  in  space  is  subdivided  into  a  set  of  smaller  sub- 
domains  or  ’’finite  elements”.  An  approximate  solution  is  then  developed  over  each  of  these. 
This  subdivision  of  a  whole  into  parts  allows  accurate  representation  of  complex  geometries  and 
nonlinear  material  behavior,  for  instance  [30]. 
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2. 3. 3.2  Choice  of  elements 

Numerous  types  of  finite  elements  are  available  within  general  purpose  finite  element  codes. 
The  choice  of  elements  for  a  particular  application  will  depend  not  only  upon  the  more  obvious 
matching  of  its  shape  attributes  to  those  of  the  physical  structure  to  be  modelled,  but  also  the 
less  obvious  matching  of  its  displacement  degrees  of  freedom  (DOF)  to  be  compatible  with  the 
structure’s  probable  displacement  behavior. 

For  example: 

Although  three-dimensional  (3D)  elements  can  be  assembled  to  model  any  structure,  care  should 
still  be  taken  to  ensure  that  several  layers  of  elements  are  present  in  regions  subject  to  bending 
moments  in  order  that  meaningful  stress  values  are  obtained  (Gebbeken  &  Wanzek  [18], 
[35]).  A  mixture  of  elements  within  a  model,  each  suited  to  a  particular  geometric  shape,  is 
preferable  to  attempting  to  use  the  same  element  types  throughout.  Thus,  for  example,  slim  two- 
dimensional  (2D)  rectangular  elements  with  four  nodes  should  not  be  used  where  slender  beams 
in  bending  would  be  preferable.  This  is  because  the  four  noded  element  possesses  only  a  linear 
displacement  interpolation  between  its  nodes.  The  element  strain,  which  is  the  differential  of  the 
displacement,  is  therefore  constant  throughout  the  element.  Thus,  if  in  reality  an  appreciable 
strain  gradient  exists  along  the  length  of  the  element,  serious  incompatibilities  will  exist.  The 
choice  of  a  higher  order,  quadratic  or  cubic  interpolation  element  (with  extra  face  nodes)  will 
considerably  reduce  this  problem  since,  for  instance,  the  quadratic  element  (with  mid  side  nodes) 
will  allow  linear  variation  of  strain  and  stress  along  and  across  an  element  to  exist. 

2. 3. 3. 3  Element  restrictions 

After  the  choice  of  suitable  elements  has  been  made  and  a  solution  has  been  obtained,  the 
analyst  should  be  aware  of  possible  deficiencies  in  the  output  information.  Table  2.2  lists  a 
range  of  elements  which  are  commonly  employed  in  FE  analysis.  The  table  lists  the  degrees 
of  freedom  or  directions  in  which  translational,  U,  and  rotational,  R,  displacements  or  forces 
are  available  for  an  element  together  with  a  list  of  the  direct,  s,  and  shear,  t,  stresses  produced 
by  the  element.  The  remarks  column  comments  upon  the  limitations,  for  example,  the  stresses 
which  are  not  given  by  the  element.  It  may  be  seen,  for  instance,  that  shell  elements  do  not 
record  the  presence  of  direct  stresses  or  shear  stresses  through  their  thickness.  This  will  be 
accurate  enough  for  thin  shells  or  for  a  thick  shell  with  a  large  radius  of  curvature.  However  it 
may  prove  to  be  too  restrictive  in  the  case  of  a  thick  shell  with  a  tight  radius  of  curvature.  In 
this  case  solid  element  modeling  might  well  be  considered. 
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Table  2.2:  Range  of  elements  commonly  employed  in  FEM  (Prinja  [?]) 


Standard  Element  Type 

Active  DOF 

Output  Stresses 

Remarks 

3D  Solid 
(8-noded  or 

20-noded  bricks.) 

Ux,Uy,Uz 

sx,  sy,  sz,  txy,  txz,  tyz 

Bending  may  require  more 

than  one  8-noded  element 

through  thickness  or  use 
higher  order  20-noded  bricks. 

Axi-sym  solid 
(3-noded  or 

4-noded) 

Ur,Uz 

sr,  sz,  shoop,  trz 

Bending  may  require  more 

than  one  3  or  4  noded  el. 
through  thickness  or  use 
higher  order  6-noded 

or  8-noded  elements. 

3D  Shell 

Ux,Uy,Uz,Rx,  Ry,Rz 

sx,  sy  and  txy  in  plane 

Check  validity  for  thin/ 

/thick  shell  assumption. 

No  direct  or  shear  stress 

through  thickness. 

Axi-sym.  Shell 

Ur,Uz,  Rrz 

smeridional  and  shoop 

No  direct  or  shear 

through  thickness. 

2D- Plane  strain 

Ux,Uy 

sx,  sy  and  txy  in  plane 

and  sz 

Through  thickness 

strain  is  zero. 

2D-Plane  stress 

Ux,Uy 

sx,  sy  and  txy  in  plane 

Through  thickness 

stress  is  zero. 

Plate 

Uz,Rx,Ry 

sx,  sy  and  txy  in  plane 

and  t  transverse 

Membrane  stress 

is  zero. 

3D-Beam 

Ux,Uy,Uz,Rx,  Ry,Rz 

saxial  ,  t  torsion 

No  transverse  shear  stress 

Truss 

Ux,Uy,Uz 

saxial 

Only  axial  stress 

2. 3. 3. 4  Stress  computation 

As  the  finite  element  method  is  an  approximation  method,  it  is  important  to  know  that  stresses 
are  dependent  on  the  location  where  they  are  calculated  (Reddy  [29]).  It  often  happens  (for 
isoparametric  elements)  that  stresses  are  most  accurate  at  Gauss  points  of  a  quadrature  point 
one  order  less  than  that  required  for  full  integration  of  the  element  matrix;  especially  for  shear 
stresses. 

For  example,  in  a  plane  stress  analysis  the  horizontal  and  vertical  displacements  are  continuous, 
and  in  the  analysis  of  a  plate  bending  problem  using  the  transverse  displacement  as  only  unknown 
variable,  this  displacement  and  its  derivatives  are  continuous,  too.  However,  this  continuity  does 
not  mean  that  the  element  stresses  are  continuous  across  element  boundaries  (Bathe  [2]).  Of 
course,  this  problem  diminishes  when  using  a  finer  mesh. 
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2. 3. 3. 5  (Further  aspects) 

Further  aspects  related  to  FE  and  error  are: 

•  Locking  effects, 

•  integration  points, 

•  calculation  of  stress  out  of  strain,  (etc.) 

(...  they  will  not  be  discussed  here  ...) 

2. 3. 3. 6  Errors  in  FE  solutions 

(taken  from  Reddy  [30]) 

The  errors  introduced  into  the  finite  element  solution  of  a  given  differential  equation  can  be 

attributed  to  three  basic  sources: 

1.  Domain  approximation  error,  which  is  due  to  the  approximation  of  the  domain.  In  two- 
or  three-dimensional  problems  involving  non-rectangular  domains,  domain  approximation 
errors  are  inevitably  introduced  into  the  FE  solution.  In  general,  these  can  be  interpreted 
as  errors  in  the  specification  of  the  data  of  the  problem  bacause  we  are  solving  the  given 
differential  equation  on  a  modified  domain.  As  we  refine  the  mesh,  the  domain  is  more 
accurately  represented,  and,  therefore,  the  boundary  approximation  errors  are  expected 
to  approach  zero. 

2.  Quadrature  and  finite  arithmetic  errors,  which  are  due  to  the  numerical  evaluation  of  in¬ 
tegrals  and  the  numerical  computation  on  a  computer.  When  finite  element  computations 
are  performed  on  a  computer  round-off  errors  of  numbers  and  errors  due  to  the  numerical 
evaluation  of  integrals  are  introduced  into  the  solution.  In  most  linear  problems  with  a  rea¬ 
sonably  small  number  of  total  degrees  of  freedom  in  the  system,  these  errors  are  expected 
to  be  small. 

3.  Approximation  error,  which  is  due  to  the  approximation  of  the  solution.  This  error  is 
inherent  to  any  approximation  method.  It  can  be  shown  that  the  approximation  error 
is  zero  for  the  single  second-order  and  forth-order  equations  with  element-wise-constant 
coefficients  ([2],  [30]). 
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2.3.  SOURCES  OF  ERRORS  WITHIN  THE  METHOD 


2.3.4  The  finite  difference  method  (FDM) 

{taken  (and  modified)  from  W.G.  Gray  [23]} 


2. 3.4.1  Introduction 


In  elementary  calculus  the  derivative  is  usually  defined  as  the  difference  between  values  of  a 
function  at  two  points  divided  by  the  distance  between  the  points  taken  in  the  limit  as  the 


distance  approaches  ”0”: 


df  f{x  +  Ax)-f{x) 

—  =  limAx-^o - - 


(2.7) 


FD  procedures  involve  the  inverse  of  the  above  process  such  that  a  derivative  based  on  a  dif¬ 
ferential  distance  dx  is  approximated(!)  as  a  difference  based  on  some  finite  distance  Ax.  Of 
course  in  making  the  approximation  a  certain  error(!)  is  introduced  (section  2. 3. 4. 2). 


FDM  is  based  on  replacing  the  differential  equation  which  holds  for  the  entire  domain  under 
study  with  a  number  of  discrete  approximations  to  the  differential  equation  which  model  the 
equation  at  selected  points.  The  mentioned  selection  of  points  within  this  discrete  approximation 
is  commonly  understood  as  discretization.  Thus  special  focus  will  be  taken  on  this  matter; 
see  later. 


Another  basic  of  this  method  is  the  dependency  of  each  node  on  it’s  neighbours.  Here  a  certain 
degree  of  dependency  on  the  nearby  points  has  to  be  taken  into  account.  A  simple  approximation 
often  used  is  that  one  point  may  only  depend  on  the  immediate  neighbouring  points.  A  complex 
approximation  conceivably  could  depend  on  all  other  points.  Thus  although  FDM  utilizes  local 
approximations,  it  is  possible  for  points  throughout  the  entire  domain  to  contribute  to  the  local 
approximation. 


A  simple  example,  one-dimensional: 

An  approximaton  on  an  Id-interval  \xi^2  <  2:  <  for  example,  is  formulated  as  follows: 

where  f  is  the  approximation  to  a  function  /.  In  this  special  example,  all  nodes  between  Xi-2  and 
Xi+i  (four  nodes)  have  been  taken  into  account.  Thus,  three  derivatives  of  f  can  be  formulated: 


df 


E 


1 


m=— 2;  m^j 


■  i^i+j  ^i+m)  f._ 


n 


-2;  k^j,m 


(x  -  Xj+k)  . 
(^i+j  ^i+k) 


(2.9) 
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1 

E 


n 


{X  -  Xi+k) 


i=^2  m=^mV:j  n=-2;  n^j,m  fc=-2; 


(2.10) 


f/3r  1  1 

=E/<+i  E 


dx^ 


m=i'm^j  n=-2^#j,m  fe=-2;^j,m,n 

/^  ^  w  \ 


E 


(2.11) 


Obviously,  the  approximations  to  the  first  and  second  derivatives  depend  on  space.  This  means 
that  different  approximations  will  be  obtained  by  different  locations.  The  third  derivative  de¬ 
pends  not  on  space  and  we  obtain  a  constant  value  when  using  a  four  node  approximation  as 
presented  here. 

Therefore,  the  accuracy  of  the  approximation  increases  by  taking  more  nodes  into  account.  And 
it  decreases  when  we  go  to  higher  derivatives  with  an  approximation  using  a  finite  number  of 
nodes.  Thus,  we  need  a  greater  number  of  nodes  for  higher  derivatives. 

When  considering  equidistant  (xj+i  —  Xj  =  Ax)  nodes  and  evaluating  the  derivatives  at,  Xi_2 
for  example,  we  obtain: 


di,  ,  -ll/(xi_2)  +  18/(xi_i)-9/(xi)-k2/(xHi) 

x,:_o)  =  - — - 


dx 


6Ax 


(2.12) 


dx^ 

dx® 


{Xi-2)  = 


{xi-2)  = 


2/(xi_2)  -  5/(xi_i)  4-  4/(xj)  -  f{xi+i) 
Ax2 

-f{xi-2)  +  3/(xi-i)  -  3/(xj)  -h  /(xj+i) 
Ax® 


(2.13) 


(2.14) 


[23] 


Numerical  examples 

We  have  selected  four  arbitrary  functions  as  numerical  examples;  fig.  2.7  plots  their  curves. 
Their  derivatives  have  been  evaluated  for  the  exact  solution  and  for  an  FD-approximation. 

Tab.  2.3  sumarizes  the  results  comparing  the  exact  solution  and  the  FD  approximation.  For  all 
functions  two  discretizations  have  been  conducted,  for  four  nodes  each: 

1. )  Ax  =  1.0;  from  Xi_2  —  2  to  Xj+i  =  5,  and 

2. )  Ax  =  2.0;  from  Xi-2  =  1  to  x^+i  =  7. 


2.3.  SOURCES  OF  ERRORS  WITHIN  THE  METHOD 


All  these  results  have  been  evaluated  at  the  first  node  Xi-2  that  is  2  for  Ax  —  1.0  and  that  is  1 
for  Ax  =  2.0. 

For  the  functions  /2  and  /a,  there  is  no  difference  between  the  exact  and  the  FD  solution.  /4 
and  fs  have  some  significant  differences  instead. 


Conclusion 

The  accuracy  (see  section  1.3.2)  of  FD  approximations  depends  on 

•  the  regularity  of  mesh, 

•  the  function  approximated, 

•  the  locations  where  the  derivatives  are  evaluated. 


•  the  distance  Ax  that  represents  the  mesh  size. 
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Table  2.3:  Selected  functions  and  their  derivatives:  Exact  solution  and  FD  approximation 


function  and  derivatives 

exact  solution 

FD  approximation 

error 

(f,  r\ry,x^xi^2 

(f,  r,  f”);  ^  =  ^i-2 

in% 

/2  :  y 

-  -(a; -4.5)2  +  18 

x—2 

Ax  —  1.0 

x=2 

-  f = 

5.0 

5.0 

0 

...  ^(*1-2)= 

-2.0 

-2.0 

0 

...  ^{xi-2)  = 

0 

0 

0 

h-  y  =  \ 

-  3x^  +  12x  —  4 

x—2  /  x—1 

x=2  /  x=l 

x=2  /  x=l 

...  £{Xi-2)  == 

2.4  /  6.6 

2.4  /  6.6 

0/0 

...  ^{xi-2)  = 

-3.6  /  -4.8 

-3.6  /  -4.8 

0/0 

...  ^{Xi-2)  = 

1.2  /  1.2 

1.2  /  1.2 

0/0 

h-- 

h-  (^3:  +  1) 

x=2  /  x=l 

x=2  /  x=l 

x=2  /  x=l 

T' 

to 

II 

4.88  /  8.56 

4.92  /  5.96 

<1%  /  30% 

...  j^{xi-2)- 

-3.12  /-4.08 

-2.64  /-4.96 

15%  /  21% 

...  ^{Xi^2)^ 

-0.48  /  -1.44 

0.96  /  0.24 

(±0 

fs: 

y  —  X  ■  sin{^)  +  9 

x=2  /  x=l 

x=2  /  x=l 

x=2  /  x=l 

•  ••  ^{^i-2)  = 

1.382  /  0.918 

1.432  /  1.395 

4%  /  52% 

0.758  /  0.120 

1.05  /  -0.038 

38%  /  (±!) 

...  ^{Xi^2)  = 

-0.469  /  -0.766 

-0.489  /  -0.626 

4%  /  18% 

•  the  number  of  nodes  that  have  been  taken  into  account, 

•  the  locations  of  these  points,  and 

•  the  order  of  derivative. 


In  contrast  to  the  FE  method  where  one  is  looking  at  a  problem  from  the  global  point  of  view, 
the  FD  method  relies  on  a  specific  number  of  local  approximations. 


2. 3. 4.2  Error  estimate 

(taken  from  Gray  [23]) 

No  error  estimate  is  obtained  directly  when  using  polynomial  interpolation  to  generate  the 
difference  expressions.  But  it  is  possible  to  determine  the  error  by  applying  the  difference 
expression  successively  to  polynomials  of  higher  degree  until  an  error  is  detected.  Then  the 
dependency  of  this  error  on  the  grid  size  and  the  derivatives  of  the  polynomial  may  be  obtained. 
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This  method  in  general  reveals  the  truncation  error  E  of  the  method.  This  can  be  well  formulated 
in  the  special  example  from  above  by  simplification  of  equations  2.12-2.14; 

(6a^-6a  +  l)^;jV  (2.15) 

12  dx^ 

Herein  a  represents  a  certain  position  between  the  interpolation  points,  where  a  can  be 

choosen  such  that  0  <  a  <  1  and  h  —  Xi+i  —  Xi. 


2. 3. 4.3  Quadrilateral  elements 

Of  course  there  are  lots  of  formulations  to  approximate  the  state  variables  with  respect  to 
elements.  In  FD,  most  commonly  a  quadrilateral  element  approach  is  used;  especially  when 
using  wave  propagation  codes:  Hydrocodes  (section  3). 

A  central  difference  scheme  ([23],  [2],  [30],  ...)  is  usually  adopted  within  the  FD  method.  The 
elements  are  commonly  defined  as  quadrilaterals.  It  is  assumed  that  the  state  variables  (p,  Sjj, 
p,  i)  are  constant  within  an  element. 


Figure  2.8:  Quadrilateral  element:  Location  of  variables 

Deformations,  velocities  and  accelerations  are  evaluated  at  each  corner  node.  Pressure,  density, 
internal  energy  and  the  stress  are  parameters  that  are  evaluated  in  the  cell  center  (fig.  2.8).  That 
means  that  the  integration  of  the  stress  divergence  over  the  element  within  the  conservation  of 
momentum  (chapter  ’Hydrocodes’  3,  equ.  3.1),  for  example,  is  simplified  by  assuming  that  the 
stress  is  constant  over  the  element.  Therefore,  very  small  element  sizes  are  demanded  when 
using  the  FD  method. 

2. 3. 4.4  Boundary  and  symmetry  problems 

Especially  the  FDM  includes  problems  at  boundaries  and  at  the  axis  of  revolution  when  using 
axial  symmetry.  This  has  been  described  in  detail  by  Benson  in  [3]  with  special  respect  to  wave 
propagation  simulations. 
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2.3.5  Comparison  of  FEM  and  FDM  -  an  overview 

Sections  2.3.3  and  2.3.4  have  demonstrated  that  finite  element  and  finite  difference  methods 
are  numerical  procedures  which  arise  from  different  governing  equations  but  which  have  more 
similarities  in  final  application  than  casual  inspection  might  reveal. 

Before  making  a  comparison  of  both  methods,  a  numerical  analyst  must  somehow  ensure  that 
he  is  using  each  method  to  best  advantage.  A  comparison  of  a  particular  FE  method  with  a 
particular  FD  method  will  not  yield  preference  of  FEM  or  FDM  in  general.  Discussions  which 
argue  that  one  method  is  better  than  the  other  in  fact  often  reflect  that  the  author  has  merely 
utilized  and  programmed  one  method  better  than  he  has  the  other  or  he  in  fact  prefers  one 
method  from  his  specific  educational  point  of  view. 

But  each  method,  FDM  and  FEM,  has  got  it’s  special  features  which  may  be  desirable  for  a 
particular  application.  Therefore,  tab.  2.4  should  provide  a  reasonable  basis  for  categorizing  a 
method  as  primarily  FD  or  FE. 

By  this  at  hand,  the  user  must  decide,  which  method  is  best  to  apply  for  his  specific  task. 
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2.3.  SOURCES  OF  ERRORS  WITHIN  THE  METHOD 


Table  2.4:  Categorizing  FEM  and  FDM  (W.G.  Gray  [23]) 


Finite  Element  Method  (FEM) 

Finite  Difference  Method  (FDM) 

Approximates  the  functions  themselves. 

Approximates  derivatives  of  functions  directly. 

Provides  global  approximations  which 
are  restricted  to  neighborhoods  of 
the  grid  to  whatever  degree  desired. 

Provides  approximations  at  points  which  borrow 
information  from  neighboring  points  to 
whatever  degree  desired. 

Structure  of  the  method  readily  lends 
itself  to  irregular  curved  grids. 

Traditionally  uses  a  regular  mesh; 
often  only  even  a  rectangular  one. 

Approximations  developed  from  families 
of  basis  functions  and  integration 
over  the  grid. 

Approximations  developed  via  curve  fitting 
for  Taylor  series  expansions. 

Resultant  difference  expression  provided 
implicitily  to  the  code  which  performs 
the  integrations.  Approximation  is 
buried  in  the  integrals. 

Resultant  difference  expression  explicitly 
provided  by  the  programmer  to  the  code. 
Approximation  my  be  ’’seen”. 

Nodes  included  in  discrete  approximations 
determined  by  basis  functions,  integration 
rule  and  weighting  factions. 

Complete  freedom  to  selected  nodes  to  be  used 
in  a  discrete  approximation  provided. 

Even  on  a  regular  grid,  different  approx, 
may  be  made  at  different  nodes 
(e.g.  9-node  Lagrangin  elements). 

On  a  regular  grid,  similar  approximations 
are  made  at  different  nodes. 

Difference  expressions  are  often  written 
in  terms  of  function  values  and  deriva¬ 
tives  of  the  functions. 

Difference  expressions  are  generally  in  terms 
of  function  values  only. 
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Table  2.5:  (continued) 


Finite  Element  Method  (FEM) 

Finite  Difference  Method  (FDM) 

Though  typically  weighted  over  different 
coordinates,  derivatives  may  sometimes  be 
lumped  by  selection  of  an  appropriate 
nodal  integration  formula  (e.g.  trapez.- 
rule  for  bilinear  elements). 

Though  not  typically  done,  derivatives  with 
respect  to  one  independent  variable  may  be 
weighted  over  another  variable  (e.g.  a  time 
derivative  approx,  may  be  a  sum  of  time  der¬ 
ivatives  evaluated  at  different  spatial  locations) . 

A  staggered  grid  is  virtually  impossible 
because  of  grid  irregularities. 

A  staggered  grid  is  easily  implemented 
because  of  regular  rectangular  nodal  placement. 

Moving  boundary  problems  are  conceptio- 
nally  straightforwarde. 

Moving  boundary  problems  are  traditionally 
difficult. 

Functional  implanting  may  minimize  need 
for  fine  mesh  near  steep  gradients. 

Large  gradients  require  a  fine  mesh  for 
resolution.  More  elements  are  required. 

Data  input  often  is  complex  and 
undiscovered  data  errors  can  be  a  cause 
of  trouble. 

Data  input  may  be  simple  because  of 
regularity. 

2.4  Simplifications  concerning  material  modeling 

Depending  on  the  real  material  behavior,  we  have  to  select  the  proper  material  model.  Material 
behavior  can  be  described  as: 

•  linear  elastic, 

•  non-linear  elastic, 

•  elasto-plastic, 

•  viscous, 

•  and  others  (and  combinations). 

Choosing  the  wrong  model  is  obviously  a  source  of  error. 

Fig.  2.9  shows  some  basic  reological  models  related  to  their  specific  material  behavior. 
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behavior 

model 

c 

1  linear  elastic 

/  Hooke  - 

^  e  -  spring  element 
(1  parameter) 

E 

a 

j 

^  ideal  plastic 

^  r 

Saint- Venant  - 

(1  parameter) 

o 

i 

f 

k  elasto-plastic 

/  y  Prandtl 

^  ^  s  fHooke  -  St.-Venant) 

(2  parameters) 

E  ^ 

^ 

a‘ 

^  viscous,  visco-elastic 

1  Newton  -  Voigt-Kelvin 

- ►  t  damping  device 

Hooke) 

(2,  3  parameters) 

E  ^ 

- -  -  ) 

Figure  2.9:  Behavior  and  models  of  basic  material  formulations 


When  integrating  material  models  into  finite  methods,  it  is  common  practice  to  simplify  in 
certain  points.  Although  simplifications  are  usual  and  always  done,  one  should  be  aware  of  the 
consequences.  Nevertheless,  there  are  unavoidable  simplifications  also  in  material  modeling. 

The  following  examples  are  not  a  complete  review  of  all  approaches  that  might  be  possible  in 
the  field  of  material  modeling,  but  they  are  implemented  in  almost  every  commercial  code  and 
are  broadly  used. 

We  have  to  insinuate  that  most  users  have  not  realized  the  consequences  resulting  from  this.  To 
understand  us  correctly,  it  is  not  a  goal  here  to  disrepute  these  methods.  But  doing  simplifica¬ 
tions  as  presented  now,  one  has  got  a  more  or  less  ”  rough  estimate”  of  the  problem  to  be  solved, 
before  any  calculation  has  been  started. 
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2.4.1  Example  1:  Hooke’s  law  of  elasticity 


Up  to  a  certain  load  level  -  the  elastic  limit  the  material  response  can  be  understood  as  ideally 
elastic.  In  this  case  Hooke’s  law  of  elasticity  is  commonly  used.  This  approach  is  necessarily 
based  on  the  homogenization  principle  (chapter  1,  fig.  1.4). 

Hooke’s  law  can  be  written  as: 

o'ij  —  Eijki  e/j/  ,  (2.16) 

wherein  Eij^i  is  the  elasticity  tensor. 

This  general  formulation  contains  81  independent  material  parameters  for  an  elastic,  anisotropic 
material. 


By  progressively  conducted  simplifications  the  number  of  81  independent  parameters  in  eq.  2.16 
can  be  reduced  to  2  which  leads  to  an  isotropic  formulation: 
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(2.17) 


The  (only)  two  paramaters  of  this  elastic  law  for  isotropic  material  are  called  Lame-constants 
A  and  /i.  They  depend  on  values  necessarily  to  be  derived  from  experiments,  the  modulus  of 
elasticity  E  and  the  shear  modulus  G;  or  the  modulus  of  compression  K  and  the  PoiSSON’s 
ratio  ly. 


E{l-v) 

(l  +  j,)(l_2zy)  ’ 

(2.18) 

7-, 

{l  +  u){l-2v)  • 

(2.19) 

Ev 

{l  +  u){l-2u)~ 

(2.20) 

^  E  Ell  —  Ei2 

2{l  +  v)~  2 

(2.21) 

The  procedure  of  simplification  to  get  this  law  is  shown  in  [17],  for  example,  tab.  2.6  gives  an 
overview  of  the  steps  that  lead  to  different  material  laws;  from  the  general  elasticity  law  to  the 
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two-parameter  model  that  is  used  by  most  numerical  codes  for  isotropic  elasticity. 


Table  2.6:  Material  definitions  and  number  of  independent  material  parameters  [17] 


material 

number  of  independent 
material  parameters 

arbitrary:  no  simplification 

81  (36) 

anisotropic 

21 

monocline 

13 

orthotropic 

9 

transversal  isotropic 

5 

isotropic 

2 

Annotation:  Most  numerical  codes,  especially  hydrocodes,  offer  only  the  simple  elasticity  law 
with  two  parameters.  As  we  know,  concrete  is  a  highly  anisotropic  material  and,  therefore,  the 
material  model  for  concrete  should  have  at  least  21  independent  material  parameters  (tab.  2.6). 
If  we  simplify  the  concrete  material  by  homogenization  techniques,  we  should  be  aware  of  the 
possible  error. 

2.4.2  Example  2:  PoiSSON’s  ratio 

Reflecting  the  Poisson’s  ratio  leed  to  a  change  in  theory;  [32]. 

The  PoiSSON-ratio  u  can  be  calculated  from  experimental  measurements  taken  from  the  uniaxial 
tension  test. 


1/  := 


ly  Afa; 


(2.22) 


-I'ex 


e.  =  -ve. 


(2.23) 

(2.24) 


wherein 

v  :  Poisson’s  ratio, 

^  :  Poisson’s  constant. 


steel  ^ 

1/ 


For  examle,  the  static  values  are: 


0.3 


(2.25) 
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^concrete  _  q  3  (0.12  -  0.46) 


(2.26) 


While  steel  is  highly  homogenious,  concrete  is  absolutely  not  homogenious,  and  the  ingredients 
change  due  to  the  ordered  concrete  quality.  If  someone  models  a  highly  homogenious  material 
the  quantity  of  the  Poisson’s  ratio  is  almost  a  constant.  Modeling  concrete,  and  setting  v  =  0.2 
means  producing  an  input  error  up  to  240%. 


Transformation  into  a  stress-explicit  description  yields: 


Inversion  yields: 


Change  in  volume: 


ai 

^ 

(2.27) 

a2 

ei,2  =  =  -2/—  , 

(2.28) 

as 

ei,3  =  =  -I/—  , 

(2.29) 

ei  =  ^[ai-  i'{a2  +  as))  , 

(2.30) 

62  =  ^[0-2  -  17{ai  -l-crs))  , 

(2.31) 

€S  =  ^[as-  v{ai  -l-cr2))  , 

(2.32) 

(2.33) 

(1 +  .)(!- 2.)  X’  >')«  +  K£i  +  «)I. 

(2.34) 

"*>=(1 +  .)(!- 2.)  XI  + 

(2.35) 

dV  -1-  iAdV  —  dx{l  +  ei)  dy{l  +  €2)  dz(l  -1-  ea)  . 

(2.36) 

Hence,  the  volume  strain: 

AdF  dV  +  MV  -dV 
®  dV  ~  dV 

_  dx{l  -I-  Cl)  dy{l  -I-  62)  dz{l  H-  £3) 
dx  dy  dz 

—  Ci  -f  €2  +  63  +  eie2  +  €2^3 

^lin&ar  ^nonlinear 

The  volume  dilatation: 


(2.37) 


- (ci  +(72  +  as)  =  — - —  Sam  +  Sara  , 


^  (^linear  — 


E 


E 


(2.38) 
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where 

o-m  =  ^(o-l +  cr2 +0-3)  (2.39) 

=  mid.  (average)  stress  =  compression  stress: 

e(=  eiinear)  —  >  (2.40) 

wherein  K  is  the  modulus  of  compression. 

For  K  follows  that  Poisson’s  ratio  only  can  reach  values  between  0  <  v  <  0.5.  Figure  2.10 
shows  the  functional  of  the  modulus  of  compression  K  over  v.  v  -A  0.5  yields  K  -A-  infinity. 


Gm  =  constant 
e  =  0 

^  K  -A  oo 
ly  -A  0.5 


Figure  2.10:  Modulus  of  compression  K  vs.  u 


Conclusion: 

For  ly  -A  0.5  theory  must  change.  For  u  is  close  to  0.5  means,  physically,  approaching  incom¬ 
pressible  materials.  This  requires  a  constitutive  law  for  the  pressure,  the  hydrostatic  part.  So 
we  must  extent  the  classical  theory  of  continua. 
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Consequence:  Split  of  the  stress  tensor 

Because  of  this,  many  methods  are  based  on  the  additive  split  of  the  stress  tensor  into  a  hy¬ 
drostatic  part  and  into  a  deviatoric  part.  The  first  describes  the  change  in  volume  and  -  as  we 
know  from  experimental  studies  (on  metals)  -  has  no  influence  on  plastic  yielding.  The  latter  de¬ 
scribes  the  shape  deformation  part.  Applying  this  for  concrete  needs  an  extention:  While  steel 
or  rubber  is  almost  incompressible,  concrete  is  compressible  due  to  the  collapse  of  concrete’s 
material  matrix  by  applying  high  pressures  (p>20  kbar).  Thus,  the  hydrostatic  part  has  to  be 
considered,  accompanied  by  an  Equation  of  State  (EoS). 


Chapter  3 


Wave  propagation  codes  - 
Hydrocodes 

3.1  Introduction 

3.1.1  Available  Hydrocodes 

Nowadays,  numerous  commercial  hydrocodes  are  available  for  almost  everyone  (tab.  3.1).  Au- 
TODYN  is  one  of  them,  and  used  by  the  authors. 

As  explained  before,  one  should  always  be  aware  of  the  method  that  is  used  by  the  software. 
Otherwise,  a  meaningful  use  is  impossible.  The  user  is  always  confronted  with  the  key  question 
’’which  software  is  the  best  for  my  problem”. 

Software  was  mostly  developed  for  special  purposes,  and,  therefore,  there  is  no  general  purpose 
software  available.  Because  of  the  competition  between  software  developers  &  vendors  they 
often  promise  too  much  (section  1.3.4). 

But,  it  is  the  user’s  duty  to  ensure  that  his  numerical  result  and  that  the  interpretation  of 
the  result  (section  3.2)  is  correct.  Not  knowing  the  ’’secrets”  of  the  source  code  worsens  this 
problem. 
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Table  3.1:  Available  hydrocodes,  taken  from  Meyers1994  [26],  modified  and  supplemented 


hydrocode 

properties 

developer 

TOUDY 

Lagrange-processor 

Sandia  Nat.  Labs. 

HEMP 

Lagrange-processor 

Lawrence  Livermore  Nat.  Labs. 

STEALTH-2d,3d 

Lagrange-processor 

Science  Appl.  Inc. 

PRONTO-2d,3d 

Lagrange-processor 

? 

MESA-2d,3d 

Euler-processor 

Los  Alamos  Nat.  Labs. 

PAGOSSA-3d 

Euler-processor 

Los  Alamos  Nat.  Labs. 

J0Y-3d 

Euler-processor 

Lawrence  Livermore  Nat.  Labs. 

DYNA-2d,3d 

Lagrange-processor 

Lawrence  Livermore  Nat.  Labs. 

LSDYNA2d,3d 

Lagrange-,  Euler-processor 

Lawrence  Livermore  Nat.  Labs. 

CALE-2d 

Lagrange-processor 

Lawrence  Livermore  Nat.  Labs. 

CAVEAT 

? 

Los  Alamos  Nat.  Labs. 

CTH-2d,3d 

Lagrange-processor 

Sandia  Nat.  Labs. 

PICES-2d,3d 

Lagrange  und  Euler  gekoppelt 

Physics  Inernational 

CRALE-2d 

Arbitrary  Lag. /Euler-processor 

? 

CSQ  II-2d 

Euler-processor 

Sandia  Nat.  Labs. 

EPIC-2d,3d 

Lagrange-processor 

Honeywell 

NIKE-2d,3d 

Lawrence  Livermore  Nat.  Labs 

ZEUS 

Lagrange-processor 

Segletis  &  Zukas 

Autodyn2d 

Lagrange  und  Euler  gekoppelt 

Century  Dynamics 

AutodynSd 

Lagrange-,  Euler-processor 

Centnry  Dynamics 

TDL  MADER 

Lagrange-processor 

C.  Mader 

DYSMASS 

Lagrange-,  Euler-processor 

IHBC 

SHARK 

Euler-processor 

3.1.2  Why  and  when  hydrocodes  are  used 

There  is  obviously  an  increasing  demand  by  building  owners  to  design  their  structures  for  hazard 
loadings.  These  loadings  are  caused  by  explosion,  impact,  penetration  and  perforation,  for 
example,  and  produce  high  pressures  in  the  kilobar  region  within  microseconds  in  the  material. 

Structural  dynamics  methods  are  well  known  but  insufficient  for  hyper  dynamics  (section  1.4). 
Thus,  wave  propagation  codes  (hydrocodes)  are  nowadays  applied. 

It  is  not  a  trivial  task  to  numerically  predict  damage  and  fracture  caused  by  hyper  dynamic 
impacts.  The  fields  of  interest  where  hydrocodes  are  applied  are  -  for  example 
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•  fluidal  flow  simulations, 

•  gas  dispersion, 

•  free  explosions  in  air, 

•  contact  explosions  on  structures, 

•  impact,  penetration  and  perforation  phenomena. 

Some  examples  will  be  shown  in  the  sections  bellow  (sections  3.5.1  and  3.5.2). 

Actually,  hydrocodes  are  the  most  powerful  numerical  tools  for  the  numerical  study  of  shock  and 
impact  phenomena.  But  up  to  now,  it  is  impossible  to  totally  replace  all  experimental  investi¬ 
gations  by  (cheaper)  hydrocode  simulations.  Experiments  are  still  necessary,  for  benchmarking, 
for  calibration  of  material  parameters,  etc.  . 

3.2  Pre-  and  post-processing,  interpretation  of  results 

In  order  to  use  finite  methods  effectively,  software  tools  must  be  user-friendly.  Therefore,  user 
graphical  interfaces,  pre-  and  post-processors  have  been  developed.  But  user-friendly  does  not 
mean  that  pre-  and  post-processors  can  replace  the  knowledge  and  experience  of  the  user. 

3.2.1  Pre-processing 
Overview: 

•  Generation  of  geometrical  models, 

•  generation  of  numerical  models, 

•  possibly  data  transfer  from  CAD, 

•  implementing  material  data, 

•  specify  units,  methods,  etc., 

•  generation  of  meshes, 

•  check  of  elements  and  input. 


•  check  of  numerical  models. 
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3.2.2  Post-processing 

The  treatment  of  numerically  derived  data,  their  sorting,  combining  and  presentation/visualization 
in  graphical  form  is  referred  to  as  post-processing.  In  the  post-processing  phase,  results  may  be 
produced  in  the  following  forms;  see  ”How  to  -  Interpret  Finite  Element  Results”  by  Baguley 
&  Hose,  NAFEMS  publications: 

•  displaced  shape  plots, 

•  contour  plots  of  displacements, 

•  plots  showing  reaction  forces  and  moments, 

•  listing  of  stresses  or  strains  sorted  by  ascending  or  descending  magnitude, 

•  contour  plots  of  stresses  or  strains, 

•  vector  plots  showing  the  directions  of  principal  stresses, 

•  plots  of  stresses  or  strains  along  paths, 

•  plots  of  deflection,  stress  or  strain  against  time, 

•  path  integrals  such  as  /-integrals  for  fracture  mechanics. 

3.2.3  Interpretation  of  results 

The  interpretation  of  results  is  a  focal  point! 

Results  from  post-processing  help  the  analyst  to  understand  how  the  model  behaves.  To  achieve 
the  objective  of  the  analysis  it  is  necessary  to  translate  these  results  into  the  behavior  of  a  real 
structure.  This  may  require  adjustments  to  the  results  to  allow  for  the  differences  between  the 
model  and  the  real  structure,  and  this  process  is  referred  to  as  interpretation. 

Increasing  computer  technology  in  connection  with  sophisticated  pre-  and  post-processors  offers 
both,  new  chances  and  new  risks. 

Besides  calculations  of  huge  problems,  the  chance  is  to  lay  a  powerful  tool  in  the  user’s  hand 
which  enables  ’’easy”  pre-processing  and  an  extensive  analysis  of  the  results. 

Risks  lie,  for  example,  in  mistakes  of  pre-  (input  problems)  and  post-processing  (misinterpreta¬ 
tions)  made  by  the  user  and/or  clients.  There  is  also  a  significant  risk  to  use  a  wrong  model  or 
the  wrong  method.  Using  tools  as  black  boxes  prevents  insight  into  what  one  is  actually  doing 
which  increases  this  problem. 
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The  code  developer  should  clearly  state  how  the  program  treats  the  data;  e.g.  ’’shaping”  or 
’’smoothing”  of  data  at  integration  points  into  contour  plots  etc.  .  Traceability  and  transparency 
(section  1.3.2)  of  tools  are  necessary  for  the  evaluation  of  methodical  inherent  errors. 

3.3  What  is  ’’code  validation”? 

3.3.1  Introduction 

Scientists,  code  users  and  developers  have  of  course  recognized  that  validation  and  verification 
is  crucial  for  the  reliability  and  safety  of  numerical  simulations.  Within  the  last  two  decades  the 
term  ’’code  validation”  has  been  used  increasingly  (see  ’’Preface”  at  page  1). 

But  code  validation  has,  however,  become  neither  a  common  term  nor  does  it  even  have  a  stan¬ 
dardized  and  commonly  understood  content.  ’’Code  validation”  is  approached  in  different  ways 
and  in  different  intensity;  depending  on  the  specific  problem  to  be  solved  (here:  hyperdynamic 
problems),  the  method  used,  the  different  goals  of  the  interested  parties  (section  1.3.4),  the 
education  of  users  and  developers,  and  the  available  budget. 

For  example.  Hills  and  TruCANO  describe  in  [24]  a  statistical  validation  of  engineering  and 
scientific  models.  Herein  the  focal  point  concerning  code  validation  is  the  match  between 
experiment  and  numeric  after  plotting  their  results.  The  key  question  that  Hills  et  al.  approach 
is:  ’’When  is  the  agreement  between  experimental  measurement  and  model  prediction  good 
enough?” 

3.3.2  Ho^v  to  find  a  proper  ’definition’ 

Validation  is  a  confirmation,  through  the  provision  of  objective  evidence,  that  the  requirements 
of  a  specified  intended  use  or  application  have  been  fulfilled  (section  1.3.2). 

Here,  the  key  requirement  is  to  fulfill  the  essential  accuracy  of  numerical  results.  Therefore,  we 
propose  the  following: 

1 .  Plausibility  studies  of  various  problems  in  order  to  know  the  limit  of  applicability  (section 
3.3.3). 

2.  Testing  codes  against  numerous  (critical)  credible(!)  experimental  results:  Benchmarking 
(section  3.3.4). 

3.  Testing  codes  against  each  other.  Do  not  use  only  one  commercial  code! 
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4.  Plausibility  testing  of  the  results  by  hand  calculations,  known  formulars  (if  it  is  possible). 

But  objective  evidence^  that  numerical  results  are  safe  and  reliable,  cannot  be  ensured  only  by 
the  user  nor  only  by  the  code  developers  each  on  their  own.  Beyond  the  proper  interaction  of  all 
interested  parties  (section  1.3.4)  it  is  also  important,  that  code  developers,  code  vendors,  code 
users  and  their  partners  validate  each  other. 

3.3.3  A  specific  plausibility  study 

The  following  might  be  done  by  the  user  as  a  first  step  of  code  validation: 

For  example,  characteristic  material  data  has  been  reliably  obtained  by  a  set  of  flyer  plate 
impacts:  EoS  :  Us  =  a  +  b-Up  (e.g.  [?]).  Thus,  a  and  b  are  the  parameters  in  the  shock  velocity 
-  particle  velocity  relation. 

The  result  of  a  numerical  simulation  -  using  exactly  the  same  conditions  (materials,  ...)  as  in 
the  experiment  -  will  also  lead  to  a  shock  velocity  -  particle  velocity  relation:  EoS*  :  Us  = 
a*  +  b*  ■  Up. 

Such  a  study  can  be  understood  as  a  combination  of  the  first  two  items  of  the  ’definition’  above 
(section  3.3.2). 

”To  all  intents  and  purposes”,  the  experimental  data  (a  and  b)  and  the  numerical  result  (a*  and 
b*)  should  not  differ  from  each  other.  But  an  exact  match  between  testing  and  numeric  {a  =  a*, 
b  =  b*)  is  only  an  ideal  case  that  we  cannot  expect.  Indeed,  a  (b)  and  a*  (6*)  will  differ  most 
likely. 

But  why? 

The  answer  to  this  question  is  the  answer  to  ”How  can  we  conduct  code  validation?’’ . 

Code  validation  is  not  ’’only  matching  EoS  with  EoS*”  and  ”be  satisfied  when  the  difference  is 
of  about  less  than  5%”  (for  example  done  by  curve  fitting) .  Proper  code  validation  is  the  answer 
to  the  question  ’’where  do  the  5%  come  from?”.  Sections  3.4  and  3.5  try  to  give  answers. 

3.3.4  An  overall  data  base  board  for  benchmarking 

The  following  proposal  for  an  overal  data  base  board  is  partly  taken  from  the  ”DFL  Data  Base 
for  RCS  Computational  Code- Validation”  [13]. 

In  our  opinion  it  could  be  well  transfered  to  our  code  validation  problem.  For  example  an  overall 
database  board  for  benchmarks. 

Thus,  it  has  been  modified  and  supplemented  here  for  the  general  purpose  to  be  discussed. 
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A  data  base  board  is  intended  to  streamline  the  process  of  establishing  confidence  in  a  given 
computational  code.  The  data  base  contains  experimental  as  well  as  numerical  results.  The  data 
has  to  be  for  many  targets,  from  very  simple  problems  to  more  complex  objects.  Complexity 
might  rise  by  both  larger  setups  and  more  complicated  problems  like  hyper  dynamic  impacts. 
Users  and  code  developers  should  be  able  to  download  the  numerical  values  of  the  data  as  needed 
establishing  the  accuracy  of  a  new  code  and  in  testing  its  limitations.  Code  users  wishing  to 
establish  confidence  in  a  given  code  can  obtain  measured  data  for  comparison  with  their  own 
discretizations  and  computations. 

We  should  create  user  groups  of  common  interest  beyond  the  limit  of  usergroups  for  a  specific 
software  package. 

Further,  users  and  developers  have  to  be  invited  to  contribute  their  computations  and  measure¬ 
ments  to  the  data  base.  The  intention  is  that  a  user  wishing  to  examine  the  evidence  supporting 
the  accuracy  of  a  certain  code  can  look  into  the  data  base  and  find  computations  with  that  code 
compared  with  measurements  for  various  targets,  including  the  details  of  the  discretization. 
This  may  save  new  users  the  expense  of  extensive  code  validation  studies  of  their  own,  as  such 
studies  are  readily  available.  Further,  a  new  user  can  examine  the  details  of  the  discretizations 
available  on  the  data  base  to  learn  to  build  his  own  models  effectively. 

Another  objective  is  to  allow  a  user  to  compare  the  performance  of  various  codes  to  choose 
the  one  most  suitable  for  a  given  application.  Thus,  if  a  given  target  has  been  solved  using  a 
number  of  different  codes,  the  user  can  examine  the  predicted  by  the  various  codes  in  relation 
to  one  another  and  to  the  measured  results.  If  such  comparisons  are  available  for  a  variety 
of  targets,  the  user  has  quick  access  to  a  comprehensive  comparison  of  the  performance  of  the 
various  computational  techniques  available  for  that  problem. 

Those  developing  a  code  may  wish  to  use  this  data  base  to  demonstrate  their  code’s  performance 
relative  to  other  codes.  By  contributing  data  to  data  base  computed  with  their  code  for  various 
problems,  the  capabilities  of  the  code  are  exhibited,  relative  to  other  methods.  Potential  users 
shopping  for  a  code  can  evaluate  various  codes  one  against  the  other.  There  is  an  advantage  for 
both  for  users  and  for  developers. 

A  contribution  to  the  data  base  has  to  have  four  key  parts: 

•  identification. 


•  problem. 
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•  method  and 

•  data. 

The  problem  being  addressed  must  be  clearly  identified,  either  as  one  already  on  the  data  base 
or  as  a  new  problem.  The  method  used  for  solving  the  problem  must  be  fully  described.  A 
certain  instance  has  to  proof  every  ingoing  contribution  before  distributing  it.  Queries  from  the 
board  chair  are  imperative  there.  There  also  should  be  a  detailed  questionnaire  to  be  filled  out 
by  one  who  wants  to  send  a  contribution.  In  this  kind  of  form  also  general  information  might  be 
asked  for;  e.g.  information  about  the  computer  resources  including  the  disc  space  needed,  the 
memory  requirements,  and  the  running  time  on  a  specific  processor.  When  measured  data  is 
submitted  to  the  data  base  a  description  of  the  measurement  facilities,  the  instrumentation  and 
the  methods  used  have  to  be  included.  To  get  information  about  experimental  results  without 
an  detailed  description  of  the  method  used  is  for  nothing!  The  problem  geometry  must  be 
described  in  detail,  including  a  drawing  if  required  for  clarity.  Measured  data  is  desirable  if  the 
problem  is  to  be  used  as  a  reference  for  code  validation,  although  computations  alone  should  be 
posted  as  an  invitation  to  others  to  solve  the  geometry  with  their  code  and  compare  with  the 
results  already  posted  for  code  validation  purposes. 

[13] 


Specific  references 

First  of  all,  the  ”DFL  Data  Base  for  RCS  Computational  Code- Validation”  of  the  Canadian  Space  Agency  has 
to  be  mentioned.  The  text  above  is  partly  taken  from  this  organisation. 

[http://lucas.dfl.crc.ca/rcs/rcs.html] 

The  ”CFD  CODE  VALIDATION  DATABASE”  provides  a  focus  for  CFD  code  validation  in  areas  deemed  most 
important  to  U.S.  aerospace  industry. 

[http:  / /iadt. arc. nasa.gov/database/ldescribe.htmlj 

’’The  development  of  standards  for  code  validation  and  verification.” 

[http:  / /www. aiaa.org/information  /  technical  /  rev-as.html] 

Institute  for  Mathematics  and  its  Applications  (IMA),  Sandia  Nat.  Lab.:  ’’Code  validation  as  a  reliability  problem” 
by  T.  G.  TVucano. 

[http://www.ima.umn.edu/reactive/abstract/trucanol.htmll 

Sandia  Nat.  Lab.:  ’’Description  of  the  Sandia  Validation  Metrics  Project”  by  Trucano  et  al.  [34] 
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3.4  The  basics  of  hydrocodes 

In  this  section,  first  of  all,  a  brief  overview  of  the  theoretical  background  of  hydrocodes  will 
be  given  (subsection  3.4.2).  Afterwards  important  demands  are  done  that  are  essential  for  all 
codes  used  (subsection  3.4.3).  Subsections  3.4.4-3.4.7  treat  our  experience  with  hydrocodes  and, 
therefore,  some  specific  proposals  are  made. 

3.4.1  Waves,  compression  waves,  shock  waves 

The  theory  of  waves  is  important  to  understand  also  hyper  dynamic  phenomena.  But  this  will 
not  be  discussed  in  detail  here.  For  such  basics  we  refer  to  [26],  for  instance. 

3.4.2  Theoretical  background 

Basically,  it  is  important  to  consider  all  aspects  mentioned  in  chapter  2  when  talking  about 
finite  methods. 

The  most  hydrocodes  are  based  on  the  FD  method  (section  2.3.4)  for  time  discretization,  espe¬ 
cially  when  shock  wave  phenomena  are  the  key  subject.  In  addition  they  are  in  general  ’’explicit 
codes”  with  respect  to  the  solution  in  time  (explicit  integration) . 

There  are,  however,  some  additional  aspects  with  hydrocodes.  While  FE  codes  are  usually  based 
on  the  direct  solution  of  the  equilibrium  equations,  hydrocodes  are  commonly  based  on  the  direct 
solution  of  the  conservation  equations  for  mass,  momentum  and  energy.  Thus,  their  governing 
equations  are  different.  The  conservation  formulas  have  to  be  fulfilled  simultaneously  in  every 
time  step  during  the  computation. 

In  order  to  well  understand  the  development  of  shock,  eq.  3.1  includes  the  equations  for  the 
hydrostatic  case.  The  pressure  is  the  predominant  action  here.  Eq.  3.1  can  also  be  applied  to 
the  one-dimensional  stress  state.  Eq.  3.1  is  valid  for  all  situations  where  the  deviatoric  part  can 
be  neglected  compared  to  the  hydrostatic  pressure  part.  This  is  common  practice  when  solid 
material  under  high  pressure  is  under  consideration.  Therefore,  eq.  3.1  has  been  derived  by 
directly  using  the  instantaneous  shock  condition  which  is  shown  in  fig.  3.1. 

Mass  p(,Us  -  Up)  =  po  Us  (a)  ^3 

Momentum  p  +  p  {Us  —  Up)^  =  Po  +  Po  Ug  (b) 

Energy  e  +  ^  (Us  -  a?)^  =  eo  +  ^  Ug  (c) 

The  initial  parameters  (po,  Po,  eo)  are  known.  The  atmospheric  overpressure  po  is  very  small 
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Ug=  shock  front  velocity 

p^=  initial  pressure 
p  =  initial  density 
e®  —  initial  spec,  enthalpy 

p  =  hydrostatic  pressure 
p  =  density 

e  =  spec,  internal  enthalpy 
Ujj=  particle  velocity 

Figure  3.1:  Shock  front 

in  comparison  with  the  maximum  pressure  and  is  usually  neglected.  The  initial  internal  energy 
iQ  —  cq  —  ^  corresponds  to  the  actual  temperature.  Finally,  the  four  unknowns  (p,  p,  e{i),  Up) 
remain  in  three  independent  equations.  That  yields  a  lack  of  one  equation.  In  general,  any 
additional  relationship  between  two  unknown  parameters  can  be  taken  into  account.  Usually  a 
formula  between  pressure,  density  and  energy  is  taken  that  is  called  Equation  of  State  (EoS) 
(e.g.  the  authors  in  [19]-[22],  [31],  [32]). 

There  are  some  problems  concerning  the  EoS.  In  this  paper,  it  is  impossible  to  go  into  detail 
here.  Nevertheless,  one  should  be  aware  that  the  EoS  can  only  be  determined  by  experimental 
investigations. 

Eqs.  3.1  (a)-(c)  reveal  that  there  is  need  to  provide  an  additional  constitutive  material  law  of 
the  form  a  —  cr(p,  e,  e).  The  material  strength  is  usually  modeled  by  a  failure  surface  /  = 
It  will  be  derived  from  the  theory  of  continuum  mechanics. 

For  further  information  an  theoretical  details  of  hydrocodes  we  refer  to  relevant  literature 
(e.g.  [1],  [3]). 

3.4.3  Important  demands  to  the  hydrocode  used 

The  following  demands  are  essential  for  all  numerical  kernels  of  codes,  but  especially  for  wave 
propagation  codes: 

•  Asymptotic  convergence, 

•  robustness, 

•  numerical  stability, 


•  accuracy. 
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•  consistence. 


3.4.3. 1  Asymptotic  convergence 
—  Discretization  errors  in  finite  methods  — 

Again,  finite  methods  are  approximate  methods  (section  2.3.2).  Thus,  methodical  inherent  errors 
have  to  be  detected  and  reduced  to  an  acceptable  minimum.  This  is  absolutely  indispensable  in 
order  to  obey  the  physical  laws.  If  one  does  not,  the  numerical  procedure  might  smear  effects 
and  afterwards  one  is  unable  to  analyse  and  evaluate  the  numerical  results  with  respect  to  the 
causes.  In  the  worst  case,  the  obtained  numerical  results  are  nonsense  and  just  numbers  free 
to  any  interpretation.  Quantifying  the  value  of  an  ’’acceptable  minimum”  is  dependent  on  the 
kind  of  problem. 

Every  finite  method  has  methodical  inherent  errors  of  discretization,  spatial  (mesh  sensitivity) 
and  in  time  (timestep).  In  order  to  produce  safe  and  reliable  numerical  results,  it  is  important 
to  determine  the  methodical  inherent  aspects.  In  section  3.5,  a  detailed  convergence  study  will 
be  discussed,  based  on  the  convergence  theorem: 


If  a  partial  differential  equation  is  discretized  by  means  of  finite  elements  or  finite  differences, 
the  numerical  solution  must  monotonically  converge  to  the  exact  solution  of  the  mechanieal 
model  if  and  only  if  the  element  size  and  the  timestep  tend  to  zero. _ 


This  convergence  theorem  is  based  on  the  assumption  that  the  selected  finite  elements  work  ... 
(no  locking  effects,  etc.). 

The  main  goal  is  to  find  a  numerical  model  which  is  as  coarse  as  possible  and  as  fine  as  necessary. 
This  can  be  a  time  consuming  effort.  Because  of  time  pressure,  if  a  comparative  study  has  been 
done,  often  just  2  elementations  are  compared.  With  two  values,  one  can  generate  a  linear 
formulation  but  not  a  monotonic  convergence.  But  there  is  an  exception:  If  the  results  of  two 
really  different  elementations  do  not  differ  more  than  3%  the  ’’rough”  mesh  must  be  seen  as 
acceptable. 

To  find  the  reasons  of  method  inherent  problems  concerning  convergence,  stability,  and  sensi¬ 
tivity  the  theoretical  background  of  the  tool  used  has  to  be  studied. 

Important  to  know  and  to  accept  is  the  fact  that,  methodical  inherent,  there  might  exist  some 
regions  where  no  convergence  can  be  adopted.  An  example  is  given  in  section  3.5. 
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3. 4. 3. 2  Robustness 

Robustness  or  smooth  behavior:  Small  variations  must  cause  small  changes  in  results!  There¬ 
fore,  sensitivity  studies  for  all  relevant  numerical  parameters  have  to  be  taken  into  account. 

The  controll  of  robustness  is  important  especially  for  hyper  dynamic  problems  where  the  vari¬ 
ation  of  the  state  variables  is  immense  in  a  short  duration  of  time.  An  example  is  shown  in 
section  3.5.2  where  internal  numerical  parameters  in  penetration  problems  are  discussed. 

In  contact  detonation  problems  this  becomes  evident  when  using  Euler-Lagrange-coupling.  Not 
only  mesh  sizes  or  the  mesh  size  ratio  (discussed  in  section  3. 5. 1.4)  influence  the  results.  Also 
the  alignment  between  the  Euler  and  the  Lagrange  mesh  is  important  to  consider.  Using  whole- 
numbered  ratios  cause  changes  in  pressure-time  plots  that  are  physically  not  explainable  ([32]). 
Therefore,  we  recommend  a  pre-study  on  robustness  for  the  special  problem  that  has  to  be 
solved.  This  has  been  done  by  the  authors  in  [32],  for  example,  for  contact  detonation  problems 
(section  3. 5. 1.4). 

3. 4. 3. 3  Numerical  stability 

Numerical  stability  (of  explicit  (hydo-)codes,  numerical  damping,  and  others).  Normally  this  can 
be  assumed  to  be  satisfied  in  commercial  codes.  Anderson  describes  the  handling  of  numerical 
stability  in  detail  in  [1] . 

3.4. 3. 4  Accuracy 

Accuracy  is  here  referred  to  the  internal  method’s  truncation  error;  see  equation  2.1,  discussed 
in  section  2.3.1.  A  sufficient  accuracy  of  commercial  hydrocodes  has  to  be  postulated.  (— >■  see 
also  ’’conformity”  in  section  1.3.2) 

3. 4.3. 5  Consistence 

Consistence  of  time  and  geometrical  discretizations  has  to  be  postulated. 

Consistence  has  to  be  understood  with  respect  to  the  choice  of  timestep  At  and  element  size 
Ax.  The  timestep  and  the  element  size  are  dependent  to  each  other.  Consistence  signifies  that 
for  At  — ^  oo  and  for  Aa;  -A'  oo  the  difference  A  can  mathematically  expressed  by  the  differential 
coefficient  d.  In  other  words  the  numerical  integraton  scheme  becomes  the  exact  solution  of  the 
differential  equation  and  the  analytical  solution  of  the  mechanical  model  will  be  approached. 
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3. 4. 3. 6  Benchmarks 

Benchmarks  are  an  important  part  of  code  validation  (section  3.3). 

By  setting  up  a  benchmark  problem  we  can  observe  that  different  users  achieve  quite  different 
results  by  using  different  methods  and  codes.  Also  different  problems  occur. 

The  benchmark  tests  should  be  conducted  numerically  and  experimentally.  But  numerical 
results  can  not  be  fitted  to  experimental  results  simply  by  non-physical  curve  fitting.  In  such  a 
case  one  cannot  learn  anything  from  the  numerical  simulation. 

It  has  become  fashionable  to  use  numerical  tools  along  with  the  experimental  investigation.  But 
sometimes,  there  exists  a  confusion  about  the  difference  between  experimentally  determined 
model  parameters  (e.g.  for  the  yield  curve  within  the  material  model)  and  tuning  the  numerical 
input  data  until  the  numerical  result  fits  the  experimental  one. 

Thus,  it  is  important  to  formulate  standard  benchmark  tests  against  which  codes  are  tested 
(section  3.3.4). 


3.4.4  How  to  determine  the  proper  timestep 


Most  hydrocodes  use  an  explicit  integration  scheme  that  leads  always  to  problems  near  bound¬ 
aries  or  the  axis  of  revolution  because  forces  and  displacements  are  usually  not  conjugated  like 
in  FE  codes.  The  calculation  advances  from  time  fi  to  time  without  any  iterations  in 
between  and  the  differences  of  the  constitutive  values  are  based  on  the  timestep  At.  Because 
shocks  travel  distances  of  millimeters  per  microsecond  (e.g.  PETN:  9000^):  the  timestep  must 
be  small  enough  to  resolve  the  temporal  details.  To  get  the  right  value  for  the  proper  timestep, 
different  criteria  have  to  be  checked.  One  of  the  most  important  is  related  to  the  smallest  mesh 
size  Ax  and  is  called  COURANT  criterion: 


At  < 


Ax 


(3.2) 


where  c  is  the  soundspeed,  which  is  in  case  of  a  shock  wave  propagation  replaced  by  the  velocity 
of  the  shock.  The  wavespeed  c  is  not  constant,  it  changes  as  the  material  changes  (density, 
pressure,  stiffness): 


c  = 


replaced  by 


(3.3) 


The  wavespeed  c  can  be  identified  in  the  p-p-plane:  It  is  proportional  to  the  variable  slope  of 
the  EoS. 
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Additionally,  the  mesh  size  Ax  changes  significantly  while  meshes  deform.  Thus,  the  CouRANT 
criterion  has  to  be  fulfilled  in  every  timestep. 

Eq.  3.2  should  ensure  that  the  shock  front  does  not  omit  a  cell.  That  has  been  stressed  by 
additionally  using  a  time  step  correction  factor,  a  so  called  "safety  factor”  (<  1),  by  that  the 
timestep  is  multiplied.  So,  no  disturbance  should  be  able  to  propagate  across  a  zone  in  a  single 
timestep.  In  case  of  shock  loads,  "disturbance”  means,  for  example,  the  shock  width  within  the 
mesh. 

Shock  development 

Fig.  3.2  shows  the  development  of  a  shock  that  arised  from  a  contact  detonation  on  a  concrete 
plate.  The  numerical  set  up  and  further  details  are  discussed  in  section  3.5.1. 


X 


Figure  3.2:  Range  of  shock  development  (zoom  of  a  crater  region) 

On  the  left  part  of  fig.  3.2,  the  contour-plot  of  pressure  is  shown  from  p  ==  0  to  p  =  Pmax  (here: 
78  kbar).  On  the  right  part,  the  distribution  of  pressure  along  the  axis  is  plotted. 

In  [4],  Benson  explains  that  the  shock  can  only  sufficiently  be  treated  if  it  is  transfered  by 
approximately  four  to  six  elements  and  the  density  of  the  mesh  must  be  fine  enough  so  that  the 
six  element  shock  width  is  small  in  comparison  with  the  dimensions  of  the  problem. 

In  the  example  presented  here  we  have  —  7%.  This  aspect  can  be  proved  by  looking 

at  the  range  of  the  shock  development  which  is  about  2.2cm  and  includes  about  six  elements 
(fig.  3.2). 

Benson’s  demand  has,  obviously,  been  satisfied  here  by  using  a  Lagrange-mesh  size  of  <  4mm 
in  the  ’critical’  zone.  Taking  into  account  that  the  grain  size  of  concrete  is  up  to  32mm,  one 
could  start  to  argue  about  the  homogenization  principle,  because  the  element  size  is  clearly  less 
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than  the  the  grain  size. 


Proposal: 

The  COURANT  criterion  has  to  be  automatically  controlled  by  the  code  used.  In  addition,  for 
the  sake  of  reliable  and  safe  numerical  simulation  of  shock  waves  with  hydrocodes,  we  propose 
a  compliance  with  Benson’s  six  element  shock  width  demand  as  a  minimum  requirement. 


3.4.5  Numerical  treatment  of  shock  oscillations  -  artificial  viscosity 


It  is  a  common  strategy  to  smear  out  the  jump  discontinuities  associated  with  shocks  by  special 
numerical  techniques.  This  is  necessary  in  order  to  properly  handle  shocks  in  a  hydrocode 
simulation.  A  standard  approach  is  the  introduction  of  an  artificial  viscosity  that  has  to  be 
incorporated  to  the  governing  equations.  The  role  of  these  terms  is  to  damp  out  the  postshock 
oscillations  (diffusion  terms).  This  is  important  to  obtain  physically  meaningful  solutions. 

There  exist  lots  of  different  expressions  for  viscosity  terms,  for  example  the  viscosity  formulation 
by  VON  Neumann  &  Richtmeyer  from  1950  ([28]).  Some  of  them  are  described  in  [4]. 


The  timestep 
grows: 


has  to  decrease  as  the  ’’diffusion”  term,  which  is  related  to 


At  < 


Ax 


A  +  +  ? 


;  A  =  2bl{Ax) 


dv 

dx 


+  62c. 


the  artificial  viscosity. 


(3.4) 


Herein  A  and  61,2  are  numerical  parameters  that  expand  the  commonly  used  artificial  viscosity 
term 

(3.5) 


dx  V  dx 


where  /i  is  the  common  coefficient  of  viscosity. 

In  the  used  code,  a  linear  q-term  that  provides  an  additional  condition  for  the  timestep 

Ax 


At  < 


2Bc 


(3.6) 


is  included.  It  contains  a  linear  artificial  viscosity  coefficient  B. 

These  tools  are  used  for  handling  shocks  by  spreading  the  shock  smoothly  over  several  mesh 
intervals.  With  this  at  hand  the  artificial  viscosity  q  has  to  be  incorporated  into  the  difference 
equations,  e.g.  for  the  pressure  p: 


d{p  +  q)/dx. 


(3.7) 
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The  presented  pressure  history  plots  (e.g.  fig.  3.11)  have  to  be  well  understood  in  that  manner. 
Thus,  it  is  evident  that  difficulties  occur  especially  at  the  high  instantaneous  pressures  in  targets 
near  a  detonation. 

Proposal: 

The  influence  of  artificial  viscosity  mechanisms  should  be  evaluated.  Small  variations  within 
those  mechanisms  must  cause  small  changes  in  results  according  to  robustness  that  has  been 
postulated  in  section  3. 4. 3. 2. 

3.4.6  The  multimaterial  Euler  processor 

Some  codes  use  Euler  processors  that  allow  more  than  one  material  in  a  single  cell.  This 
multimaterial  procedure  includes  certain  rules  for  mixing  the  parameters  with  respect  the  surface 
(or  volume)  fractions  of  the  different  materials.  But  that  unavoidably  leads  to  difficulties  because 
some  parameters  (e.g.  the  density)  are  mixed  depending  on  the  fractions  of  the  materials  in  the 
Euler-cell,  whereas  other  parameters  (e.g.  the  velocity)  have  to  be  treated  independently  from 
the  materials. 

This  implies  that  slip  is  impossible  on  material  interfaces  within  a  multimaterial  element.  There 
are  no  slide  lines  in  these  elements.  In  our  special  problem  (the  simulation  of  contact  explosions) 
the  flow  of  detonation-gases  next  to  the  solid  body  has  to  be  modelled.  A  boundary  layer  with 
a  nonzero  shear  stress  occurs  in  the  solid  and  the  velocity  profile  of  the  gas  is  altered  because 
of  the  ’’artificial  sticking  condition”  within  mixed  elements  (Benson  [4]). 

Proposal: 

One  should  be  aware  of  problems  concerning  the  multimaterial  processor.  Corrective  measures 
are  mesh  refinements  which  might  reduce  these  problems. 

3.4.7  Treating  distorted  cells 

The  mesh  distortion  in  the  Lagrange-mesh  during  a  simulation  of  a  contact  detonation  is  shown 
in  fig.  3.3,  where  the  region  nearby  the  crater  has  been  zoomed  out. 

It  is  evident  that  numerical  results  might  not  be  trustworthy  if  large  distortions  occur.  In 
numerical  simulations  one  has  to  take  care  that  the  quadrilateral  elements  cannot  be  distorted 
into  ’’boomerangs”  or  ’’bowties”,  as  shown  in  fig.  3.4. 


o  M  rr-iTTrr’  r>  A  or/^o  r\T?  ■LT\/r\T:>r\r^r\r\T?Q 


Figure  3.3:  Distortion  of  Lagrange-elements  nearby  the  crater 


If  the  timestep  is  small  enough  and  the  strains  are  not  too  large,  it  is  impossible  to  push  a  node 
of  a  cell  into  its  opposite  edge  because  the  numerical  volume  then  approaches  zero,  leading  to  an 
infinite  residual  pressure.  The  numerical  integration  scheme  is  only  exact  on  the  master  element 
which  is  a  square.  The  more  distorted,  the  poorer  the  result. 
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A  big  problem  of  self-  or  adaptive  meshing  is  that  elements  might  be  already  distorted  in  the 
initial  unloaded  situation.  We  don’t  use  self-  or  adaptive  meshing  here. 


Figure  3.4:  Large  mesh-distortion  and  hourglassing 


The  distortion  of  elements  can  be  controlled  by  different  strategies: 

1.  The  inner  angle  Oi  between  edges  has  to  be  <  180°. 

2.  The  determinant  of  the  Jacobian  Matrix  is  not  allowed  to  become  negative,  otherwise 
dV  —  det{J)  dr  ds  dt  yields  negative  numerical  volume. 

3.  The  change  in  ratio  of  diagonals  (1:1  for  the  quadratic  master  element),  see  fig.  3.3 

Using  ’’constant  stress  quadrilateral  elements”,  might  cause  zero  energy  modes.  This  is  called 
hourglassing  (fig.  3.4).  To  ensure  that  no  hourglassing  occurs,  certain  control- algorithms  have 
been  included  in  the  most  commercial  codes. 

About  damping  algorithms 

Handling  hourglassing  by  special  tools  is  one  specific  example  for  a  damping  algorithm.  There 
might  also  be  others,  like  smoothing  peak  pressures.  Such  algorithms  shuold  be  treated  very 
carefully.  For  a  detailed  insight  we  refer  to  Benson  [4]. 

3.4.8  Conclusion 

Section  3.4  has  given  insight  into  the  basics  of  hydrocodes.  Some  typical  problems  have  been 
revealed,  and  demands  and  proposals  have  been  stated,  for  the  purposes  of  a  quality  improvement 
of  hydrocode  simulations. 
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Nevertheless,  the  generally  stated  demands  and  proposals  have  to  be  well  targeted  to  the  special 
problem  that  has  to  be  solved  by  the  user. 

To  illustrate  and  deepen  the  background  given  in  this  section,  two  specific  hyper  dynamic  prob¬ 
lems  are  exemplarily  selected  and  discussed  in  the  following  reliability  study  (section  3.5): 

•  Contact  detonation  (section  3.5.1), 

•  penetration  and  perforation  (section  3.5.2). 
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3.5  Reliability  study 


The  following  examples  are  own  numerical  simulations  with  Autodyn  concerning  reliability 
and  safety.  Special  focus  is  done  on  asymptotic  convergence  that  has  been  well  defined  in  section 
3.4.3.I. 

3.5.1  Example  1:  Contact  detonation 

3. 5. 1.1  Load  level  regions 


Explosive  charges  cause  various  regions  of  stress  and  failure  situations,  as  shown  in  fig.  3.5. 


Figure  3.5:  Contact  detonation  on  a  concrete  plate,  regions  of  different  stress  states  and  failure  situations 

As  a  result  of  the  detonation,  regional  limited  material  failure  can  be  observed.  At  the  top 
surface  of  the  plate,  the  concrete  material  mainly  fails  due  to  high  compression.  This  crushing 
causes  a  crater.  At  the  opposite  side,  the  shock  wave  is  reflected  and  converted  into  a  tension 
wave.  Because  of  the  low  resistance  of  concrete  to  tension,  one  can  observe  a  tensile  spalling. 
Between  these  regions  there  is  a  part  where  concrete  suffers  certain  levels  of  damage. 

3. 5. 1.2  Numerical  set  up 

Fig.  3.6  shows  the  numerical  set  up  for  the  Hydrocode  Autodyn2d.  Because  of  the  problem 
we  can  take  advantage  of  the  axial  symmetry.  The  explosive  charge  of  650g  TNT  (numerically 
described  by  the  well-known  JONES- Wilkins-Lee  equation)  is  centered  on  a  circular- plate  made 
of  unreinforced  concrete.  The  charge  is  in  contact  with  the  concrete  surface. 

The  Eulerian  mesh  wherein  air  and  expanding  detonation  gases  (material  without  strength)  are 
simulated  is  coupled  with  the  Lagrangian  mesh  that  enables  the  modeling  of  the  solid  (material 
with  strength)  structure.  This  requires  an  overlapping  of  these  two  meshes.  The  choosen 
overlapping  (here  10cm)  has  to  ensure  that  the  crater-region  of  the  deforming  Lagrange-mesh 
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Figure  3.6:  Numerical  simulation  of  a  contact  explosion  on  a  concrete  plate 


is  sufficiently  overlapped  by  the  Euler-mesh  at  any  time.  Therefore,  the  proper  overlapping 
measure  can  only  be  found  with  the  help  of  a  pre-study. 

In  order  to  prove  asymptotic  convergence  (section  3. 4. 3.1)  and  to  quantify  remaining  errors  the 
problem  set  up,  which  is  given  in  fig.  3.5  and  3.6,  is  examined. 


3. 5. 1.3  How  to  choose  characteristic  parameters 

Applying  displacement  formulated  finite  elements,  displacements  converge  even  better  than  their 
derivatives,  e.g.  strains.  Because  the  stresses  are  directly  related  to  the  strains  by  constitutive 
laws  for  the  material,  the  convergence  of  stresses  or  strains,  respectively,  must  be  proved.  In 
this  case  we  can  state: 

Asymptotic  convergence  of  the  displacements  is  only  a  necessary  condition,  the  convergence  of 
strains  (stresses)  is  a  necessary  and  sufficient  condition. 

Or  more  general:  Asymptotic  convergence  has  to  he  carried  out  with  respect  to  the  parameters 
of  the  governing  equations  (characteristic  parameters),  that  include  the  highest  order  of  the 
derivative  of  the  given  problem. 

The  parameters  given  in  the  governing  equations  for  hydrocodes  (eq.  3.1)  are  pressure,  volumetric 
strain  (or  density)  and  internal  energy  (or  enthalpy)  ([1],  [3],  [4]).  These  are  parameters  on  the 
lowest  level  of  the  formulated  constitutive  equations  if  macromechanics  is  adopted.  So,  these 
parameters  are  characteristic  parameters  to  be  chosen  for  convergence  studies.  Consequently, 
the  proof  of  their  convergence  is  necessary  and  sufficient.  The  proof  of  convergence  of  integrated 
values,  e.g.  impulse,  is  only  a  necessary  condition. 

Doing  so,  the  convergence  studies  here  are  carried  out  with  respect  to  the  parameter  ’’pressure 
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p”  (figs.  3.8  and  3.10)  which  is  calculated  from  stresses. 

3. 5. 1.4  Euler-Lagrange  mesh-ratio 

The  Euler-Lagrange  coupling  algorithm  provides  a  simultaneous  use  of  both  processors  in  one 
numerical  simulation  (e.g.  [3],  [12]).  The  proper  interaction  of  the  Euler  mesh  and  the  Lagrange 
mesh  depends  -  among  other  things  -  on  the  Euler-Lagrange  mesh-ratio. 

A  detailed  mesh-ratio  study  has  been  done  by  the  authors  in  [20].  The  results  are  resumed  in 
this  section. 

The  columns  in  fig.  3.8  represent  the  maximum  pressure  ’P’  versus  the  5z,ap. /^Ewi.-ratio,  where 
Slag,  is  the  Lagrangian  and  is  the  Eulerian  mesh  size.  Ratio  2,  for  example,  means  that 
there  are  twice  as  much  Euler- elements  than  Langrange-elements  with  respect  to  the  edges 
(fig.  3.7). 

^Lagr.  _  2 
^Euler  1 


Figure  3.7:  Example:  Lagrange-Euler-Ratio  =  2:1 


Figure  3.8:  Maximum  pressure  vs.  mesh  sizes  and  -ratio  in  target  #0  (figs.  3.6  and  3.9) 

One  result  from  the  study  displayed  in  fig.  3.8  is  that  the  ratio  between  the  both  mesh  types 
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From  fig.  3.8  one  can  see  that  -  with  respect  to  the  suggested  restriction  (eq.  3.8)  -  the  quality 
of  the  result  is  almost  independent  of  the  mesh-ratio  if  eq.  3.8  is  fulfilled.  Fig.  3.8  also  reveals 
an  increase  of  peak  pressure  when  refining  the  mesh.  This  is  in  contradiction  to  the  convergence 
theorem  (section  3.4.3.1),  and  -  obviously  -  in  fig.  3.8  asymptotic  convergence  cannot  be  observed. 
This  will  be  explained  in  section  3.5. 1.5. 


Because  of  the  tremendeous  rise  of  calculation  time  with  increasing  number  of  cells,  a  ratio  of  2 
(as  shown  in  fig.  3.7)  has  been  choosen  for  our  elementation  (pre-)study.  Conducting  the  same 
study  with  a  ratio  of  3  would  multiply  the  calculation  time  with  a  factor  of  >  3,  whereas  the 
significance  of  the  result  concernig  elementation  would  be  the  same. 


Proposal: 

For  the  simulation  of  contact  detonation  problems  using  Euler-Lagrange  coupling,  a  ratio  be¬ 
tween  2  and  4  is  best  with  respect  to  both,  proper  exchange  of  data  (pressure,  particle  velocity, 
displacements  between  the  two  different  mesh  types)  and  lowest  expense  of  calculation  time. 

According  to  the  robustness  criterion  (section  3. 4. 3. 2),  we  propose  a  ratio  of  ^  =  ^±5, 

where  d  should  be  about  10%  of  the  Euler  element-size  [32]. 


3. 5. 1.5  How  to  Choose  Target  Points 

Based  on  the  physics,  the  structure  which  is  numerically  simulated  should  properly  respond 
everywhere.  Due  to  the  convergence  theorem  (section  3. 4. 3.1)  one  would  expect  convergence 
everywhere,  in  general.  That  this  might  not  be  the  case  is  illustrated  in  this  section.  The  study 
is  based  on  the  problem  given  in  figs.  3.5  and  3.6.  We  also  refer  to  what  has  been  mentioned  in 
section  3. 5. 1.4. 

Fig.  3.9  highlights  the  target  points,  already  given  in  fig.  3.6. 

We  recall  fig.  3.8  which  data  correspond  to  target  #0,  that  it  is  directly  nearby  the  contact  det¬ 
onation.  Here  the  methodical  inherent  problem  in  determining  the  discretization  error  becomes 
clear:  No  asymptotic  convergence  can  be  observed  in  target  #0. 

In  fig.  3.10  the  results  of  three  target  points  have  been  compared,  with  respect  to  the  fineness  of 
the  mesh.  We  could  not  reach  asymptotic  convergence  in  targets  #0  and  #1.  Only  beginning 
with  target#  2  asymptotic  convergence  could  be  proved.  Several  studies  with  hydrocodes  led 
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Figure  3.9:  Target  points 


to  the  answer  that  one  cannot  expect  convergence  at  targets  nearby  the  range  of  instantaneous 
shock  development.  The  reasons  are  described  in  section  3.4  and  that  is  not  specific  for  the  used 
code  only. 

Convergence  can  be  expected  if  the  target  points  are  out  of  a  nearby  region  around  the  contact 
charge.  Very  close  to  the  explosive  charge  convergence  is  not  guaranteed. 

Fig.  3.11  shows  the  time  plot  of  the  pressures  with  respect  to  different  meshes.  The  pressure-time 
plots  correspond  to  target  ^1.  One  can  see  the  different  peaks  and  shapes  of  the  pressure-time 
plots  with  respect  to  different  mesh  sizes.  The  mesh  size  varies  here  from  3.0  to  0.5  mm  in  the 
Euler  grid.  The  Lagrangian  grid  is  twice  the  size,  as  mentioned  before  (it  runs  from  6.0  to  1.0 
mm). 


3. 5. 1.6  Resume 

The  greater  the  distance  of  the  target  points  from  the  contact  charge,  the  better  the  asymptotic 
convergence  of  pressure-time  plots.  But  in  order  to  sufficiently  compute  the  stress  state  in  the 
material,  it  is  necessary  to  ensure  convergence  in  all  possible  regions  even  as  close  as  possible  to 
the  contact  charge.  If  convergence  has  not  been  proved,  one  does  not  know  how  far  away  the 
numerical  solution  is  from  the  exact  solution. 

In  case  of  nearby-detonations,  when  there  is  no  contact  between  charge  and  structure,  hydrocodes 
reproduce  the  experimental  facts  qualitatively  and  quantitatively  in  all  regions  of  the  structure, 
but  not  in  case  of  contact  detonations.  Therefore,  the  target  points  should  be  located  where  the 
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Figure  3.10:  Peak  Pressure  vs.  mesh  fineness  for  different  targets 


elements  remain  almost  undistorted,  where  the  strain  rates  are  moderate  (<  10  j)  and  where 
the  shock  wave  velocity  is  in  the  range  of  the  soundspeed. 

Of  course,  the  pressure  peaks  and  the  stress  peaks,  respectively,  are  the  most  sensitive  values, 
whereas  the  impulse  converges  faster,  because  it  is  an  integrated  value.  But  it  is  a  wrong  con¬ 
clusion  from  this  to  choose  the  impulse  as  characteristic  parameter  for  asymptotic  convergence 
studies.  The  impulse  is  not  a  parameter  in  the  constitutive  equations. 

It  would  be  very  useful  to  rezone  and/or  to  refine  the  Lagrangian  mesh  if  it  becomes  significantly 


CHAPTER  3.  WAVE  PROPAGATION  CODES  -  HYDROCODES 


74 


PRESSURE  (Mbar)  Target#  1 


Figure  3.11:  Pressure  vs.  time  for  different  mesh  sizes  for  target#  1 

distorted.  But,  the  user  has  to  pay  by  increasing  computation  time.  In  addition,  it  is  not  an 
easy  task  to  obtain  a  clear  declaration  about  the  convergence  from  the  momentum-time  plot. 
The  time-integration  of  pressure  is  unfortunately  not  unique.  The  integration  limits  are  not 
easily  choosen  because  of  the  fact  that  the  pressure  intersects  the  time  axis  somewhere,  but  not 
at  the  same  time  for  every  mesh  size.  To  choose  a  specific  time  limit  for  the  integration  of  the 
pressure,  produces  no  declaration  about  convergence  as  well,  because  the  momentum  would  not 
be  unique  by  this. 

The  momentum  that  is  derived  by  the  integration  of  the  pressure  with  respect  to  time  is  shown 
in  fig.  3.12. 

The  only  possible  statement  is  that  at  critical  target  points  (#0  and  #1)  the  momentum  does 
not  asymptotically  converge  to  a  certain  reference  momentum-time  curve. 

Conclusion: 

The  conclusion  taken  from  these  studies  is  that  one  can  not  expect  convergence  in  a  critical 
region  around  the  contact  detonation.  In  the  presented  case  of  a  specific  contact  detonation,  the 
boundary  of  this  region  is  of  about  3  cm  (8  elements)  around  the  contact  surface.  Inside  of  this 
zone  the  numerical  results  are  not  reliable  (fig.  3.13). 

But  where  are  the  reasons  for  this  ? 

A  pre-study  only  with  air  that  is  modelled  by  an  ideal  gas  Equation  of  State  (without  any 


Figure  3.12:  Momentum  vs.  time  for  different  mesh  sizes 


Figure  3.13:  The  ’’critical”  region  around  the  contact  detonation 


additional  constitutive  laws)  has  been  carried  out.  In  a  limited  region  within  an  Euler-mesh,  the 
initial  specific  internal  energy  has  been  modified  so  that  initial  pressures  similar  to  an  explosion 
have  been  approached.  Doning  so,  neither  contact  nor  coupling  of  Euler-  and  Lagrange-grids 
was  necessary.  Even  in  that  simple  case,  convergence  problems  concerning  pressure  nearby  the 
high-pressure-zone  have  occured  [32]. 

But  at  our  special  research  subject  of  contact  explosions  on  concrete  additional  problems  occur. 
So,  one  should  have  deep  knowledge  and  a  right  feel  for  all  the  problems  that  occur  especially 
at  the  overlapping  zone.  The  problems  are: 

•  the  coupling  of  Euler  and  Lagrange  (overlapping)  especially  the  exchange  of  informations. 
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•  the  nonlinearity  in  constitutive  laws, 

•  large  Lagrange-element  distortions, 

•  the  nonlinearity  of  deformation, 

•  extremely  high  pressures  and  rates, 

•  a  constant  stress  assumption  similar  to  a  single  point  integration. 

If  the  interaction  between  load  and  structure  has  to  be  considered,  the  coupling  of  Euler  and 
Lagrange  meshes  is  enforced.  If  in  addition  the  load  is  shock-like,  problems  in  the  overlapping 
zone  occur.  In  section  3.4  some  methodical  inherent  problems  are  briefly  discussed  in  order  to 
well  understand  the  reasons. 


3.5.2  Example  2:  Penetration  and  perforation 
3. 5. 2.1  On  the  numerical  set  up 


A  concrete  plate  is  impacted  by  a  steel-copper-projectile.  The  numerical  set  up  is  shown  in 
figure  3.14. 


^^^^j~^xis  of  revolution 
steel-copper-projectijle  (lagrange  mesh) 
concrete  plate  (R=20cm/T=6cm) 


lagrange  mesh 


a 

o 
■  CO 


scaled  mesh  resolution  e.g.  1mm  x  1mm  (constant) 
K - >< — ^ ^  ’ 


R=20cm 


Figure  3.14:  Numerical  setup  of  a  penetration-perforation  study 

Here  we  use  the  hydrocode  Autodyn2d  and  the  developed  concrete  model  by  Ruppert  ([32], 
[22]).  The  projectile  has  got  a  steel  kernel  and  a  copper  liner.  The  tip  has  got  an  ogive  kind  of 
shape.  For  details  we  refer  to  [32]. 
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Focus  on  the  projectile 


The  tip  has  been  meshed  by  joining  the  nose  part  and  the  body  part  (fig.  3.15).  The  given 
projectile  dimensions  are  8  •  26  [mm]. 


Figure  3.15:  Projectile  elementation 

Figure  3.15  shows  the  elementation  of  the  projectile.  One  specific  cell  of  the  nose  has  been 
zoomed  out.  This  element  is  inadequate  and  might  cause  numerical  difficulties  due  to  the 
lagrange  mesh  deformations  (section  3.4.7). 

We  run  the  risk  that  this  element  might  be  distorted  into  a  boomerang  or  into  a  bowtie.  Thus, 
also  an  alternative  meshing  method  for  the  projectile  nose  has  been  adopted  as  shown  in  fig.  3.16. 


Conclusion 

The  evaluation  of  exit  velocities  of  the  projectile  for  the  two  different  projectile  meshes  is  shown 
in  figure  3.17.  Version  2  has  got  no  ’’problem  cells”. 

The  results  (concerning  the  exit  velocity)  are  similar  but  differences  are  observable.  But  the 
results  differ  relatively  small  and  they  are  obviously  based  on  the  difference  of  mesh  fineness 
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Figure  3.17:  V-t  and  s-t-diagram 

(mesh  resolution)  between  the  two  elementations:  Version  1  has  got  16  elements;  version  2  has 
got  12  elements. 

Nevertheless,  we  propose  to  use  version  2  for  further  penetration/perforation  studies,  because 
mesh  distortion  problems  cannot  occur  at  least  at  the  beginning  of  the  calculation. 

In  the  following  section  the  inflence  of  the  mesh  resolution  for  the  concrete  plate  will  be  shown 
briefly. 

3. 5. 2. 2  Convergence 

As  mentioned  above  in  detail  (section  3.4.3. 1),  the  elementation  (mesh  density,  ...)  might  sig¬ 
nificantly  influence  the  numerical  results.  To  draw  conclusions  also  in  the  field  of  penetration 
and  perforation,  a  mesh  density  study  has  been  conducted  (see  Ruppert  [32])  and  the  results 
are  shown  here. 

It  is  important  to  determine  the  appropriate  mesh  densities  and  to  find  an  optimum  mesh  ratio 
if  possible,  in  order  to  fulfill  the  convergence  theorem  (section  3.4.3. 1). 

For  this  convergence  study,  a  concrete  plate  with  T  =  6cm  of  thickness  has  been  choosen.  The 
projectile  (metal  inside  and  copper  shell)  has  got  an  initial  velocity  of  vq  =  455^  •  Experimental 
investigations  have  not  been  carried  out  here. 

The  results  are  evaluated  by  space-time-diagrams  with  respect  to  the  exit  velocity,  that  increases 
with  mesh  density  (fig.  3.18). 

As  shown  in  figure  3.18  asymptotic  convergence  can  be  achieved:  The  exit  velocity  is  not  sig¬ 
nificantly  increasing  by  refining  the  (concrete-)mesh  from  1mm  to  0.5mm  and  again  to  0.25mm. 
A  further  mesh  refining  will  not  significantly  change  Vexit- 
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lagrange  mesh,  concrete  [mm] 


Figure  3.18:  Asymptotic  convergence  of  the  exit  velocity  for  increasing  mesh  density,  for  Vq  — 

455f 

3. 5. 2. 3  Influences  of  internal  parameters 

According  to  the  algorithm  applied,  internal  parameters  govern  the  numerical  calculation.  Pen¬ 
etration/perforation  problems  are  usually  solved  by  a  Lagrange- Lagrange  coupling  algorithm 
(e.g.  Benson  [3]).  Thus,  we  focus  on  internal  Lagrange-Lagrange  interaction  parameters  here. 
A  summary  of  our  penetration/perforation  study  is  briefly  shown  in  this  section  and  the  most 
important  internal  parameters  have  been  evaluated  (taken  and  modified  from  Ruppert  [32]). 

The  relevant  interaction  parameters  that  have  to  be  specified,  are: 

•  Gap  size  and  gap  algorithm:  Numerical  methods  based  on  FE  discretization  don’t  work 
when  cells  are  overlapping  or  penetrate  into  each  other  (section  3.4.7).  In  consequence  of 
the  Lagrange-Lagrange  coupling  algorithm,  the  involved  Lagrange  meshes  must  have  a 
certain  distance  to  each  other  in  order  to  avoid  cell  overlapping.  Therefore,  a  gap  size  has 
to  be  well  specified.  This  gap  is  not  physically,  it  is  just  a  numerical  feature!  It  can  be 
modelled  geometrically  (visible  gap)  or  it  can  be  represented  invisible  (internal  gap)  [12]. 

•  Friction:  Friction  between  materials  at  the  slide-line  of  their  meshes  can  be  considered  by 
a  specific  coefficient  of  friciton.  Friction  might  also  be  neglected  in  certain  circumstances 
depending  on  the  problem  to  be  solved. 

•  Time  step:  How  to  determine  the  proper  time  step  is  crucial  for  explicit  codes  (section 
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3.4.4).  The  time  step  is  multiplied  by  a  specific  safety  factor  (<  1).  This  prevents  over¬ 
lapping  of  opposite  cells,  because  by  this  feature  it  is  impossible  to  overjump  the  gap  in 
one  single  time  step. 

•  Erosion:  Erosion  is  a  special  feature  to  transform  elements  into  free  mass  points  by  an 
user  specified  erosion  criterion.  This  technique  might  be  powerful  to  handle  distorted  cells 
that  are  caused  by  deep  penetrations  and  perforations  [12],  For  example:  If  a  specific 
effective  plastic  strain  is  exceeded,  the  relevant  cell  will  be  eroded.  In  our  problem 
this  signifies  that  the  protectile’s  outer  layers  are  stripped  off  during  the  penetration  into 
a  concrete  plate. 

•  Inertia:  The  user  has  got  the  possibility  to  swich  on  or  to  swich  off  mass  inertia  effects 
of  the  eroded  cells  (if  erosion  is  used).  We  observed  that  the  erosion  criterion  is  mesh 
sensitive.  Thus,  erosoin  and  inertia  of  eroded  cells  have  to  be  used  carefully  [32]  and  we 
prefere  a  detaild  pre-study  when  using  these  numerical  features. 

The  variation  of  internal  parameters  from  the  default  values  by  the  user  might  cause  changes 
of  the  numerical  results.  It  is  the  user’s  duty  to  validate  their  influences.  Therefore,  primarily 
the  robustness  has  to  be  checked  (section  3.4. 3. 2). 

Our  Lagrange/Lagrange-inter action  study  reveals  the  influences  of  the  most  important  internal 
parameters,  qualitatively  as  well  as  quantitatively.  The  results  are  shown  in  tab.  3.2  where  the 
projectile’s  exit  velocity  is  the  observed  parameter. 


Table  3.2:  Influences  of  internal  Lagrange/Lagrange  parameters 


Internal  inter- 

’’default” 

Variation  done 

this  influences 

action  parameter 

(Autodyn2d) 

^exit 

Gap  size  [mm] 

0,20 

(0,08)  0,13  -  0,25 

(30%)  4% 

Gap  algorithm 

visible 

internal 

25% 

Friction 

off 

on:  0,02  -  0,08 

observable,  but  remote 

Time  step  (safety  factor) 

(3=0,2 

0,1  -  0,2 

observable,  but  remote 

Erosion  (projectile) 

off 

27% 

Inertia  of  eroded  cells 

off 

on 

25(!)% 

It  becomes  evident,  that  specific  internal  parameters  influence  the  results.  For  example,  the 
inclusion/exclusion  of  the  inertia  of  eroded  cells  or  the  value  of  the  numerical  gap  size  are 
significantly  changing  the  projectile’s  exit  velocity.  Up  to  here,  none  material  parameter  has 
been  discussed. 
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Nevertheless,  all  the  variations  done  here  are  within  the  limits  suggested  by  the  tool’s 
manuals! 

Conclusion: 

Without  this  study  no  substantial  statement  can  be  made  upon  the  choice  of  internal  interaction 
parameters  within  the  Lagrange/Lagrange- interaction  algorithm.  Thus,  we  recommend  not  to 
just  trust  on  the  defaults. 

And  the  question  comes  up:  ’’Are  all  internal  parameters  accessible  and  changeable  by  the 
user?”  Of  course  this  is  just  a  rhetorical  question  and  some  of  the  interaction  parameters 
are  hidden  behind  the  ’’black  box”.  The  influences  of  all  internal  parameters  will  not 
appear  to  the  user  unless  the  code  developers  let  him  know  the  ’’secrets”.  But  for  the  purpose 
of  quality  assurance  of  numerical  simulations,  we  attach  importance  to  the  traceability  of  all 
internal  parameters  used  in  the  codes. 

3. 5. 2. 4  Remarks 

Although  this  study  on  Vexit  is  important  to  do,  it  is  not  a  complete  and  sufficient  convergence 
study  for  penetration  and  perforation  phenomena  at  all.  The  asymptotic  convergence  of  the  exit 
velocity  is  only  a  necessary  but  not  a  sufficient  condition. 

Thus  strictly  speaking,  a  sufficient  fulfillment  of  convergence  implies  the  asymptotic  convergence 
of  the  methodical  inherent  state  variables  within  the  conservation  equations. 

Hyperdynamic  problems  solved  by  the  authors 

In  [32]  some  numerical  studies  are  given  with  the  following  topics: 

•  Free  and  near  field  detonation, 

•  indoor  explosions  +interaction  with  structure, 

•  contact  detonation, 

•  bodies’  impact,  penetration  and 


•  perforation 


Chapter  4 


Proposal  to  develop  user  guidelines 


One  never  can  stop  learning,  if  one  wants  to  keep  up  with  the  time.  This  law  becomes  obvious 
particularly  with  regard  to  using  numerical  methods  and  tools.  This  chapter  focuses  on  the  user. 
Every  user  is  livelong  a  trainee;  the  assistant  or  person  in  charge  who  executes  an  application 
and  the  scientific  user  who  develops  the  numerical  tool  as  well  as  the  person  who  must  be  able 
to  assess  and  evaluate  the  results. 

4.1  Intention 

The  intention  is  to  accompany  the  user  for  the  sake  of  a  credible  use  of  finite  methods  by 
user  guidelines.  The  theoretical  background  has  been  supplied  before  (chapter  3).  This  chapter 
outlines  the  basics  to  develop  individual  and  managable  user  guidelines  for  specific  hyperdynamic 
problems. 

4.2  The  cycle  in  general 

In  the  following  focal  points  that  have  been  mentioned  above  and  that  are  relevant  for  the  user 
are  sumarized  here: 

•  The  task 

•  Description  of  reality 

•  Mechanical  model 

•  The  method 
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•  The  tool 

•  Handling  results 

Our  conception  has  to  be  understood  as  a  checklist  that  accompanies  the  user  from  the  first 
step,  understandig  the  task  to  the  final  step:  handling  the  results. 

We  use  the  term  ’’cycle”,  because  this  schedule  has  not  to  be  treated  just  by  the  top-down 
method  beginning  from  the  first  topic  (’’The  task”)  and  ending  at  the  final  topic  (’’Handling 
results”).  Rather,  the  user  might  continually  jump  between  the  topics  back  and  forth.  For 
example,  the  topic  ’’The  tool”  can  not  handled  seperately  from  the  method  on  wich  the  tool  is 
based  on.  Also  the  topic  ’’Mechanical  model”  can  not  be  ticked  off  when  results  of  the  numerical 
simulation  are  discussed;  changing  the  mechanical  model  and  afterwards  starting  the  calculation 
again  might  often  be  advantageous.  Scrutinizing  the  results  will  often  force  the  user  to  go  the 
upper  topics. 

Of  course,  the  items  have  to  be  individually  arranged  and  adjusted  to  the  specific  task  that  has 
to  be  solved  by  the  user.  Thus,  this  list  is  certainly  not  comleted.  Depending  on  the  task,  some 
items  listet  in  the  following  are  unnecessary.  Therefore,  it’s  on  the  user  to  make  his  individual 
checklist  before  starting  to  solve  his  specific  problem. 
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4.3  Checklist 


Task 

•  Get  deep  insight:  Comprehension  of  the  problem  to  he  solved 

•  What  is  the  goal  of  the  study? 

•  List  all  interested  parties  and  their  relations  (section  1.3.4,  fig.  1.1) 

•  Classification  of  the  task  (section  1.4.1): 

—  Loading  situations  and  loading  rates, 

—  Structural  behavior  /  local  behavior  (section  1.4.3) 

—  Total  distortions 
—  Peak  pressures  (estimation) 

•  Survey  and  research: 

—  Have  similar  tasks  already  been  solved? 

—  Are  test  results  available? 

•  Evaluation/estimation  of  the  task: 

—  Application  of  resources:  Personnel  (user(s)),  tools  (software)  •  time 
—  Is  additionally  external  expertise  advantageous? 

Description  of  reality 

•  Geometry 

•  Motions 

•  Materials 

•  Loading  situations 

•  Boundary  conditions 

•  Symmetries 

•  Extent  (spatial  and  temporal) 

•  Other  environmental  conditions 

Mechanical  model 

•  Plot  all  specific  key  parameters  in  a  list: 

—  Geometry:  Extentions,  dimension  (2d/3d),  symmetry,  spatial  restrictions 
—  Structural  components:  Rods,  plates,  springs,  damping  devices,  connections 

—  Boundaries:  Constraints  and  supports,  environmental  conditions,  flow  out  of  materi¬ 
als 


85 


4.3.  CHECKLIST 


—  Loadings:  Applying,  temporal  process,  loading  combinations 
—  Materials:  Constitutive  equations,  specific  material  laws 

-  Interactions:  Describe  the  contact  of  components  (friction,  rigid  interfaces,  ...) 

•  Simplifications: 

—  Description:  List  all  simplifications  done 
—  Discussion:  Evaluate  the  influences  of  simplifications 

•  Consider  a  couple  of  possible  mechanical  models,  not  only  one 

•  see  also  section  2.2 

Method 

•  What  methods  are  alternatively  available? 

-  Element  Methods: 

*  FDM:  Finite  Difference  Method 

*  FEM:  Finite  Element  Method 

*  BEM:  Boundary  Element  Method 
—  Meshless  Methods: 

*  SPH:  Smoothed  Particle  Hydrodynamics 

*  EFG:  Element- Free  Galerkin 

*  (and  others) 

•  What  ist  the  best  method?  (section  2.3.5) 

•  Time  discretization: 

—  Explicit  or  implicit  integration? 

—  Time  step  pre- restriction:  min.,  max.,  -fsafety  factor 

-  Time  step  regulation  in  the  explicit  integration  scheme 

•  see  also  section  2.3 

Numerical  modeling 

•  Transfer  the  parameters  of  the  mechanical  model  into  the  code  (section  2.2): 

-  Are  further  restrictions  and  simplifications  necessary? 

-  Which  additional  (numerical)  parameters  are  introduced? 

*  Choice  of  a  consistent  unit-system 

*  (...) 

-  Transfer  to  the  software;  see  ’’Tool”:  Preprocessing  algorithm 

•  Spatial  discretization: 

-  Choice  of  elements:  Triangle,  quadrilateral,  element  functions 

-  Definition  of  wrapup  criteria 
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—  Are  further  pre-restriction  necessary? 

—  Are  boundary  conditions  changing  during  the  simulation? 

•  Specific  items: 

—  Internal  truncation  error  and  numerical  stability  (sections  2. 3.4. 2,  3.4. 3. 3,  3. 4. 3.4) 

—  Handling  shock  oscillations:  Artificial  viscosity  (section  3.4.5) 

—  Treating  the  multimaterial  (Euler-) processor  (section  3.4.6) 

—  Euler-Lagrange  coupling,  Lagrange-Lagrange  interaction  (sections  3.5.1,  3.5.2) 

—  Treatment  of  distorted  elements:  Handling  hourglassing,  avoid  bow  tie  (section  3.4.7) 
—  Return  algorithm:  According  to  associated  or  nonassociated  flow 
—  Damping  algorithms 
—  Fluid  transport  algorithm 

—  Stress  computation:  e.g.  a  constant  stress  asumption  within  a  single  cell  (section 
2.3.4.3) 

•  Definition  of  wrapup  criteria: 

—  Temporal:  Time  step 
—  Spatial:  Smallest  element  size 
—  Internal:  Int.  energy  vs.  reference  energy 
—  Stress  (pressure) 

—  Cutoff’s: 

*  Velocity 

*  Boundary  conditions:  e.g.  flow  criteria 

*  and  possibly  others 

•  Contact  algorithms: 

—  Interaction  procedure:  e.g.  penalty  algorithms,  gap,  exchange  of  state  variables 
—  Formulation  of  friction:  rate  dependency? 

•  Point  out  all  sources  of  errors;  see  section  2.3. 

•  (see  also  chapter  3) 


Tool 

•  What  tools  are  available?  (section  3.1.1) 

•  Limits  of  applicability  and  accordance  with  the  methods? 

•  Are  independent  calculations  possible/necessary? 

•  —>■  choice  of  the  best  tool  according  to  the  best  method. 

•  What  is  provided  by  the  tool  and 

•  assumed  to  be  axiomatic  right?  (for  example  the  JWL-EoS  for  explosives!) 

•  Input  check: 
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-  (The  handling  the  tool  according  to  the  preprocessor  should  be  ’’ticked  off’) 

-  Take  care  on  consistent  unit  systems  and  on 

-  the  material-fill  method 

•  Integrity  check: 

—  Are  all  items  available  within  the  tool,  that  have  been  mentioned  at  the  topic  The 
method”? 

*  If  so: 

•  How  are  this  items  transfered  to  and  integrated  into  the  tool  used? 

•  Are  the  concerning  parameters  changeable/variable  by  the  user? 

*  If  not  so,  initiate  further  steps: 

•  Go  through  the  manuals  in  detail,  research  and/ or 

•  ask  the  software  developers  for  more  information 

•  First  step;  approach  the  problem  by  small  pre-calculations: 

—  Define  small  and  handable  tasks  as  extracts  of  the  whole  problem 

-  (•••) 

•  Further  step;  benchmarking  (sections  3.3,  3. 4. 3. 6): 

—  Benchmarking  against  test  results 

-  Benchmarking  against  different  tools 

•  Next  step;  controlling: 

-  to  what  extent  can  the  results  be  transfered  from  the  first  step  and  from  benchmarking 
to  the  specific  problem? 

—  Region  of  divergence?  Is  the  convergence  theorem  satisfied?  (sections  3. 4. 3.1,  3.5) 

—  Conduct  the  ’’first  step  of  code  validation”  (section  3.3.3) 

-  check  Benson’s  demand  (section  3.4.4) 

-  modify  all  specific  parameters: 

*  validate  (weight)  the  parameters 

*  check  robustness  of  all  governing  parameters 

*  check  robustness  of  all  internal  parameters 

*  (for  robustness  see  also  section  3.4.3) 

•  Final  step;  simulation  of  the  task: 

-  Not  only  with  one  single  tool 

-  Not  only  by  one  single  (by  the  same)  person 
—  Not  only  one  single  mechanical  model 

—  Conduct  a  simplificated  2d-simulation  and  compare  to  3d 

•  Modificate  and  start  again  from  ’’Mechanical  model” 
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Handling  results 

•  Post- processing  and  visualization  of  results  (section  3.2): 

—  Possible  plots  are: 

*  displaced  shape  plots, 

*  contour  plots  of  displacements, 

*  plots  showing  reaction  forces  and  moments, 

*  listing  of  stresses  or  strains  sorted  by  ascending  or  descending  magnitude, 

*  contour  plots  of  stresses  or  strains, 

*  vector  plots  showing  the  directions  of  principal  stresses, 

*  plots  of  stresses  or  strains  along  paths, 

*  plots  of  deflection,  stress  or  strain  against  time, 

*  path  integrals  such  as  J-integrals  for  fracture  mechanics. 

*  Multiplots:  Mixed  plots  of  all  above 

—  Take  care  of  the  following: 

*  Plot  results  not  only  at  the  end  of  a  calculation 

*  Where  is  the  exact  position  of  the  target  points? 

*  Are  the  target  points  fixed  or  (Lagrange-)cell  related  moving? 

*  Switch  on/off  the  posprocessor’s  smooth  algorithms 

•  Comparison  to  analytical  results  (if  possible) 

•  Comparison  to  experimental  results  (see  benchmarking) 


89 


4.4.  INTEGRATION  INTO  A  QUALITY  MANAGEMENT  SYSTEM 


4.4  Integration  into  a  quality  management  system 

According  to  the  checklist  above  (section  4.3),  the  user  has  to  design  a  suited  and  individual 
checklist  concerning  his  special  problem  to  be  solved.  A  critical  evaluation  of  all  the  listed  items 
has  to  point  out  crucial  points.  This  individual  (modified)  checklist  leads  to  a  specific  user 
guideline. 

It  can  also  be  understood  as  a  quality  manual  document,  that  describes  the  total  process  in 
detail  and  states  results  achieved,  or  that  provides  evidence  of  activities  performed. 

For  the  sake  of  an  objective  evidence  (section  1.3.2),  this  quality  manual  might  be  accredited 
by  an  accreditation  council  (section  1.3.3).  The  main  characteristics  of  this  council  have  to 
be  independence,  objectivity  and  high  quality.  The  agency  has  to  be  oriented  towards  the 
realization  of  both,  academic  principles  and  practical  usage,  e.g.  for  a  proper  code  validaton  and 
benchmarking  (section  3.3).  It  has  to  be  neither  influenced  by  any  states,  governments,  interest 
groups  nor  by  lobbyists.  Like  a  scientific  institute  at  the  university,  academic  freedom  and 
autonomy  have  to  be  respected:  They  may  regulate  their  own  quality  and  standards,  but  at  the 
same  time  they  must  guarantee  transparency  of  process  and  public  accountability  in  discharging 
this  self-regulation. 

The  central  goal  is  to  provide  transparency  as  well  as  to  meet  all  the  requirements.  By  this 
at  hand,  maintenance  and  enhancement  of  quality  of  numerical  results  is  supported.  A  parallel 
objective  is  to  provide  guidance  and  information  for  beginners  (e.g.  students),  code  users,  code 
developers  and  also  scientific  institutions. 
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during  the  performance  of  this  contr^.(’ 

DATE: 

Name  and  Title  of  Authorized  Offil 
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